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Abstract 

Model  Predictive  Control  (MPC)  is  the  class  of  control  methods  that  optimizes  a 
specified  performance  index  over  a  set  of  future  inputs  to  minimize  future  output 
deviations  from  a  specified  trajectory,  subject  to  system  constraints.  MPC  operates  on  a 
receding  horizon,  calculating  a  series  of  future  control  inputs  at  each  time  step.  The 
controller  then  implements  the  first  input,  discards  the  rest,  and  calculates  the  next  series 
of  inputs  at  the  next  time  step.  Because  this  control  method  involves  on-line 
optimization,  it  has  traditionally  been  applied  mainly  to  low-bandwidth  processes. 

This  research  effort  applies  an  MPC  strategy  to  a  high-performance  aerospace 
system  with  the  goal  of  exploiting  the  constraint  handling  qualities  of  MPC  for  fault- 
tolerant  control.  To  demonstrate  the  effectiveness  of  MPC  at  handling  control  surface 
failures,  one  or  more  control  surfaces  on  a  high-performance  fighter  aircraft,  primarily  the 
F-18  High  Alpha  Research  Vehicle,  are  rendered  inoperable  during  longitudinal 
maneuvers.  In  general,  simulated  failures  of  a  single  control  element  occur  at  peak 
deflection,  whereas  the  time  and  magnitude  of  multiple  control  element  failures  are 
selected  on  a  case-by-case  basis.  Subsequent  to  the  simulated  failure  or  failures,  the 
controller  is  forced  to  compensate  by  using  the  remaining  control  surfaces  both  to 
maintain  stability  and  to  attempt  to  retain  near-nominal  performance. 

To  accomplish  these  objectives,  Matlab  and  proprietary  predictive  control  Matlab 
routines  are  used  to  develop  the  functions  and  matrices  necessary  to  implement  an  MPC 
controller.  Simulations  of  aircraft  performance  subsequent  to  control  surface  failures  are 
then  accomplished  using  Simulink  in  order  to  determine  controller  effectiveness. 
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MODEL  PREDICTIVE  CONTROL  OF  AEROSPACE  SYSTEMS 


1.0  Introduction 


1.1  Model  Predictive  Control 

Model  Predictive  Control  (MPC)  is  the  class  of  control  methods  that  optimizes  a 
specified  performance  index  over  a  set  of  future  input  moves  in  order  to  minimize  the 
weighted  future  output  deviations  from  a  setpoint  trajectory  along  a  prediction  horizon 
and  the  weighted  control  increment  inputs  implemented  along  a  control  horizon.  At  each 
discrete  time  step,  this  optimization  process  determines  a  series  of  future  control  input 
moves,  implements  the  first  move,  and  discards  the  rest.  In  the  case  of  Stable 
Generalized  Predictive  Control  (SGPC)  or  other  formulations  employing  a  stabilizing 
inner  feedback  loop,  the  output  of  the  optimization  is  a  series  of  “quasi -reference”  signals 
instead  of  control  increment  inputs  because  the  optimizer  no  longer  feeds  directly  into  the 
plant  [2,3].  This  introduction  of  the  “quasi-reference”  signal  as  the  input  to  the  inner  loop 
also  allows  for  an  independent  optimization  horizon,  r,  because  an  optimal  series  of  plant 
inputs  are  no  longer  being  directly  calculated.  In  addition,  the  inner  feedback  loop 
enforces  an  implicit  terminal  constraint,  allowing  stability  to  be  achieved  through  the 
existence  of  a  monotonically  decreasing  cost  function. 

MPC  controllers  have  traditionally  been  used  in  slow  industrial  processes  because 
they  involve  the  use  of  on-line  optimization.  However,  as  computers  become  more 


powerful,  MPC  will  consequently  become  more  practical  in  high  bandwidth  applications 
such  as  aerospace  systems.  Currently,  MPC  has  a  distinct  advantage  over  other  types  of 
control  methodologies  in  that  it  has  the  ability  to  take  into  account  hard  constraints  on  the 
system  such  as  actuator  position  limits,  actuator  rate  limits,  and  state  or  output  limits  such 
as  the  maximum  possible  normal  acceleration  of  an  aircraft.  As  mentioned  earlier,  the 
state  space  formulation  of  MPC  employed  in  this  research  effort  also  has  the  useful 
feature  of  achieving  guaranteed  stability  for  all  feasible  optimizations  based  on  the 
existence  of  a  monotonically  decreasing  cost  function. 

1.2  Importance  of  Research 

Because  MPC  is  capable  of  handling  system  constraints  such  as  actuator  rate  and 
position  limits  explicitly,  it  is  of  interest  to  employ  it  to  study  the  problem  of  aircraft 
actuator  failures.  In  the  case  of  a  damaged  or  failed  control  surface,  quick  action  by  the 
pilot  may  be  required  to  maintain  positive  control  of  the  aircraft.  However,  especially  in 
the  case  of  a  statically  unstable  aircraft  such  as  a  modem  fighter,  certain  circumstances 
may  arise  in  which  the  pilot  is  incapable  of  controlling  the  aircraft  subsequent  to  a  control 
surface  failure.  At  this  point,  the  controller  must  either  partially  or  totally  intercede  to 
prevent  a  catastrophic  failure  and  loss  of  the  aircraft. 

Given  sensor  inputs  from  the  control  surfaces  to  determine  the  existence  of  a 
failure,  an  MPC  controller  will  automatically  redistribute  control  power  in  the  attempt  to 
both  maintain  stability  and  track  the  given  setpoint.  It  should  be  noted  that  multiple 
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control  surface  failures  or  severe  out-of-trim  failures  of  certain  control  surfaces  may 
prevent  stabilization  of  the  aircraft.  Nevertheless,  depending  upon  which  control  surface 
fails  and  the  severity  of  the  failure,  an  MPC  controller  can  not  only  maintain  stability,  but 
also  maintain  near-nominal  system  performance  as  well. 

The  advantages  of  using  a  fault-tolerant  controller  are  thus  elear.  In  the  case  of  a 
severe  out-of-trim  failure  of  a  critical  control  element  (e.g.  stabilator,  thrust  vectoring 
nozzle,  etc.),  such  a  controller  will  often  be  able  to  stabilize  the  aircraft,  leading  to  the 
preservation  of  both  the  aircraft  and  the  pilot.  For  minor  failures,  a  fault-tolerant 
controller  should  be  able  to  maintain  near-nominal  performance. 

1.3  Research  Objectives 

The  primary  objective  of  this  thesis  is  to  explore  the  fault  tolerance  capabilities  of 
a  state  space  formulation  of  MPC.  Subsequent  to  a  simulated  control  surface  failure  on 
the  F-1 8  High  Alpha  Research  Vehicle  (HARV),  the  MPC  controller  will  attempt  to 
maintain  stability  and  nominal  performance  while  tracking  a  setpoint.  Secondary 
objectives  include  studying  the  effect  of  data  sample  rate  and  MPC  horizon  lengths  on 
system  performance.  These  objectives  will  be  accomplished  using  Matlab  to  develop  the 
MPC  controller  and  Simulink  to  demonstrate  its  effeetiveness  through  simulation. 

1.4  Thesis  Overview 

Chapter  2  addresses  the  essential  elements  of  developing  a  state  space  formulation 
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of  MPC.  Coprime  factorization,  plant  modification  to  accept  incremental  control  inputs, 
a  stabilizing  inner  feedback  loop,  prediction  and  constraint  equations,  and  basic  stability 
criteria  for  an  MPC  controller  are  discussed. 

Chapter  3  contains  information  about  the  F-1 8  HARV  and  the  methodology  used 
in  this  research  effort. 

Chapter  4  presents  the  results  of  the  Matlab  simulations. 

Chapter  5  offers  conclusions  based  on  the  results  of  the  simulations  and  ideas  for 
further  avenues  of  research  in  the  area  of  fault-tolerant  control  using  MPC. 
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2.0  Review  of  Literature 


2.1  Historical  Development  of  MPC 

This  section  presents  the  historical  foundations  of  the  state  space  MPC  scheme 
employed  in  this  research  effort.  Generalized  Predictive  Control  and  Stable  Generalized 
Predictive  Control  are  briefly  discussed,  with  emphasis  placed  on  their  respective 
advantages  and  disadvantages. 

2.1.1  Generalized  Predictive  Control 

The  fundamental  concept  of  Generalized  Predictive  Control  (GPC)  [1]  is  to  use  an 
explicit  plant  model  to  predict  the  plant  outputs,  y,  along  a  prediction  horizon,  N,  based 
on  inputs  implemented  over  a  control  horizon,  and  then  minimize  the  expected  sum  of 
the  Euclidean  norm  of  plant  output  deviations  from  a  setpoint  trajectory,  s,  and  the 
weighted  Euclidean  norm  of  control  activity  by  finding  an  optimal  set  of  incremental 
control  inputs.  Aw.  The  Single-Input  Single-Output  (SISO)  expectation  cost  function 
representing  this  process  is  given  by 


where  X  is  a  weight  applied  to  the  control  usage.  An  expectation  operator  is  used  because 
the  plant  model  contains  a  noise  term  colored  by  a  user-specified  polynomial. 
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GPC  has  several  advantages  over  other  forms  of  controllers.  It  does  not  require 
prior  knowledge  of  the  closed  loop  poles,  nor  does  it  suffer  ill  effects  from  singularities 
due  to  closely-spaced  poles  and  zeros  [2] .  Additionally,  GPC  allows  for  the  incorporation 
of  both  input  and  output  constraints  in  the  control  algorithm  and  provides  for  a  wide 
variety  of  tuning  parameters  that  are  useful  in  satisfying  the  given  performance 
objectives.  Unfortunately,  traditional  GPC  has  the  deficiency  of  having  neither  a  general 
guaranteed  stability  result  nor  a  general  guaranteed  robust  stability  result. 

2.1.2  Stable  Generalized  Predictive  Control 

Stable  Generalized  Predictive  Control  (SGPC)  retains  the  advantages  of  GPC,  but 
adds  guaranteed  nominal  stability  through  the  creation  of  a  stabilizing  inner  feedback 
loop  based  on  the  Youla  parameterization  of  all  stabilizing  controllers  [2].  The  closed- 
loop  poles  are  chosen  such  that  Finite  Impulse  Response  (FIR)  behavior  is  achieved.  FIR 
behavior  implies  that  the  system  will  achieve  a  steady  state  value  over  a  finite  horizon 
and  is  useful  in  that,  when  used  in  conjunction  with  a  cost  function  that  minimizes 
tracking  error,  it  can  guarantee  stability  and  asymptotic  tracking  through  the  existence  of 
a  monotonically  decreasing  cost  function  given  specific  conditions  on  the  horizons  are 
satisfied  [3].  The  MPC  optimization  is  then  wrapped  around  the  inner  stabilizing 
feedback  loop,  and  because  the  output  of  the  MPC  optimizer  no  longer  feeds  directly  into 
the  plant,  the  cost  function  must  now  be  minimized  with  respect  to  a  “quasi-reference” 
signal  to  be  described  in  Section  2.2.2. 
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2.2  Mathematical  Preliminaries  for  State  Space  MPC 

In  this  section,  the  mathematical  background  information  necessary  to  develop  an 
MPC  controller  is  presented.  This  information  includes  right  and  left  coprime 
factorization,  stabilizing  inner  feedback  loops  based  on  the  Youla  parameterization  of  all 
stabilizing  controllers,  and  plant  modification  to  accept  control  increments. 

2.2.1  Coprime  Factorization 

Although  the  development  in  [2]  is  in  terms  of  a  SISO  plant,  the  theory  can  be 
expanded  to  address  the  case  of  a  Multi-Input  Multi-Output  (MIMO)  system,  as  is  shown 
in  [4].  Consider  the  case  of  a  discrete  time  plant  Gp{z)  with  ^  inputs  and  r|  outputs.  For 
any  proper  transfer  function  matrix  with  both  a  stabilizable  and  a  detectable  realization,  it 
is  always  possible  to  describe  the  plant  in  terms  of  a  doubly  coprime  factorization  [5] 

G/z)  -N^{z)Mp\z) 

=  M-^\z)N^{z) 

where  M^{z),  N^iz),  M^iz),  iV^(z)  e  and  M^(z),  M^(z)  are  square  and  non-singular. 

For  the  purpose  of  introducing  integral  action  into  the  system,  it  is  convenient  to  modify 
the  plant  so  that  it  operates  on  incremental  changes  in  the  control  signal  instead  of  the 
absolute  control  signal.  To  do  so,  provided  that  Gp{z)  does  not  have  any  zeros  at  z  =  +1, 
the  differencing  operator  A  =  1  -  z  '  is  employed,  and  a  modified  plant  G(z)  is  defined: 
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G(z)  =  ~G^(z) 

=  N{z)M-\z) 
=  M'\z)N{z) 


Additionally,  we  shall  assign  M{z)  and  N{z)  to  be 

M(z)  =  AM^(z) 
N{z)  =  TV  (z) 


(2.4) 


where  Tl/(z),  TV(z),  M(z),  TV(z)e  ^  and  both  pairs  are  coprime. 

From  the  definition  of  coprimeness  [5],  two  matrices  M{z),  N{z)  e  are  right 

coprime  if  they  have  the  same  number  of  columns  and  there  exist  two  matrices 
X(z),  f(z)  e  having  the  relation 


X  y] 


N 

M 


XN  +  YM  =  I 


(2.5) 


Two  matrices  M(z),  N{z)  e  are  left  coprime  if  they  have  the  same  number  of  rows 
and  there  exist  two  matrices  X{z),  Y(z)  e  <R3{„  having  the  relation 


X 

Y 


NX  ^  MY  -  I 


(2.6) 
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For  any  proper  transfer  function,  it  is  possible  to  find  both  the  left  and  right  coprime 
factorizations  simultaneously  through  the  solution  of  the  generalized  Bezout  identity: 


1 

'  N  Y 

_MN 

-MX 

(2.7) 


In  a  state  space  context,  the  right  coprime  factorization  is  achieved  through  the 
introduction  of  a  state  feedback  variable,  v.  Consider  a  strictly  proper  discrete  time  plant 
G(z)  that  acts  on  incremental  changes  in  the  control  input,  ^u{k). 


G{z)  = 


A 

B 

C 

0 

(2.8) 


with  {A,B)  stabilizable  and  {C,A)  detectable.  The  state  space  representation  is 

x{k  +  1)  =  Ax{K)  +  BAu{k) 
y(k)  -  Cx(k) 

Now,  introduce  the  state  feedback  variable,  v(A:),  of  the  form 


(2.9) 


v(t)  =  z;'[  Ai/(t)  -  ] 


(2.10) 


where  is  chosen  such  that  the  eigenvalues  of  ^  +  BF^  all  lie  within  the  unit  disk.  By 
the  assumption  of  stabilizability,  this  can  always  be  accomplished.  is  chosen  to  satisfy 
a  relation  defined  in  Section  2.3.2.  Substituting  equation  (2.10)  into  equation  (2.9)  yields 
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(2.11) 


x(A:  +  1)  =  (A  +BF^)x{k)  +  BZ^v(k) 

Au{k)  =  F^x(k)  +  Zv{k) 
y{k)  -  Cx(k) 

The  transfer  matrix  from  v(k)  to  Au(k)  is  thus  seen  to  be  represented  by 


M(z)  - 


and  the  transfer  matrix  from  v(k)  to  y(k)  is 


N{z)  = 


{A  +  BF) 

BZ^ 

F 

Z 

r 

r 

ik)  is 

(A  +  BF^) 

BZ 

r 

C 

0 

(2.12) 


(2.13) 


Therefore,  the  relations  between  the  state  feedback  variable  and  the  input  and  output  are 

Au{z)  =  M(z)v(z) 

y(z)  =  N(z)v{z)  (2.14) 

>^(z)  -  N{z)  M-\z)Au{z) 


Using  the  dual  procedure  to  the  one  outlined  previously,  the  left  coprime  factorization  of 
G(z)  is  given  by: 


_M  n]: 


A  +LC 

L  B 

z,c 

Z,  0 

(2.15) 


with  a  choice  of  L  such  that  all  the  eigenvalues  of  ^  +  LC  lie  within  the  unit  disk. 
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Figure  2.1  Stabilizing  Inner  Feedback  Loop 


2.2.2  Stabilizing  Inner  Feedback  Loop 

Given  the  assumption  that  G(z)  can  be  described  by  a  coprime  factorization,  the 
transfer  function  matrices  X{z),  Y{z\  X{z),  Y{z)  e  exist  and  have  the  relation 

shown  in  equation  (2.7).  Solving  the  Bezout  identity  for  a  coprime  factorization  as  in 
equation  (2.7)  then  leads  to  the  set  of  all  stabilizing  rational  feedback  controllers  for  Giz), 
which  are  parameterized  by 

K  =  -{Y -QNYH^^QM)  (2.16) 

for  any  Q  6  [5].  Figure  2.1  depicts  the  setup  of  the  stabilizing  inner  feedback 

controller  for  a  MIMO  system,  where  the  “quasi-reference”  signal  v{k)  e  is  the  output 
of  the  MFC  optimization  routine,  y{k)  e  is  the  system  output,  and  Lu{k)  e  is  the 
incremental  control  input  such  that  Au(k)  =  u(k)  -  u(k-l).  The  particular  construction  of 
this  feedback  loop  leads  to  the  recovery  of  the  relations  in  equation  (2.14),  the  proof  of 
which  may  be  found  in  [4]  and  will  not  be  repeated  here.  One  will  note  that  the  transfer 
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function  matrix  Q  does  not  appear  in  the  recovery  of  equation  (2.14)  and  may  therefore  be 
assumed  to  be  zero.  Although  this  feedback  controller  is  implementable  with  any  choice 
of  M(z)  and  N{z),  choosing  them  to  be  FIR  operators  leads  to  a  FIR  relationship  between 
y(z)  and  Au(z)  and  the  “quasi-reference”  signal,  v(z). 

2.2.3  Plant  Modification  to  Accept  Control  Increments 

Consider  the  following  discrete  time  state  space  realization: 

xik  +  1)  =  Ax{k)  +  B^u{k)  +  Bjd{k)  ^ 

y{k)  =  Cx{k)  +  w{k-\)  +  n{k) 

where  x{k)  e  is  the  state  vector,  uQc)  e  is  the  input  vector,  d(k)  e  R'’  is  a  state 
disturbance  that  may  include  both  measured  and  unmeasured  terms,  y{k)  6  R’^  is  the 
measured  output,  w{k)  e  R^  is  an  output  disturbance  vector,  and  n{k)  represents  a 
corruption  of  the  output  vector  due  to  noise.  In  this  form,  however,  the  system  is  not 
useful  for  MPC  in  our  context  here  because  it  operates  on  the  absolute  control  signal, 
u{k),  and  not  the  incremental  change  in  the  control  signal,  Au{k).  To  remedy  this 
situation,  we  may  augment  the  plant  with  a  differencer,  which  is  the  equivalent  of  adding 
^  integrators  at  the  plant  input.  This  method  is  not  useful  in  the  case  of  non-square 
systems,  however,  and  will  therefore  not  be  utilized  in  this  research  effort.  A  more 
versatile  method  of  modifying  the  plant,  which  is  equivalent  to  adding  a  bank  of 
integrators  at  either  the  plant  input  or  output,  may  be  expressed  as 
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(2.18) 


Ax(k  +  1) 

ySk^\) 


=  A 


y(k)  =  C 


Ax(k) 

Ax{k) 

y.ik) 


+  B^Au{k)  +  B^Ad{k)  +  B^Aw{k) 


n{k) 


where 


[  ^  ol 

B 

u 

B, 

A  = 

>  - 

CB^_ 

CB, 

K  - 


0 

I 


,  C  - [ 0  /] 


(2.19) 


and  where  Ax(k)  =  x(k)  -  x(k-l),  Au{k)  =  u{k)  -  u{k-\),  Ad{k)  =  d{k)  -  d{k-\), 

Aw{k)  =  w(k)  -  ■H’(^-l)  are  the  change  in  the  state  at  time  k,  the  change  in  the  control  input 
at  time  k,  the  change  in  the  state  disturbance  vector  at  time  k,  and  the  change  in  the  output 
disturbance  at  time  k,  respectively,  is  the  vector  of  outputs  uncorrupted  by 
measurement  noise. 


2.3  State  Space  Formulation  of  Constrained  MPC 

This  section  presents  the  development  of  a  state  space  formulation  of  constrained 
MPC.  First,  the  MPC  horizons  will  be  defined,  and  their  relation  to  internal  stability  for  a 
system  operating  with  output  feedback  will  be  introduced.  Next,  the  method  for 
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developing  an  internally  stabilizing  feedback  loop  that  produces  a  FIR  relationship 
between  the  output  of  the  MFC  optimization  and  the  plant  input  and  output  will  be 
presented.  Then,  state  and  input  prediction,  input  and  output  constraint  formulation  for  a 
quadratic  optimization,  and  system  stability  will  be  discussed.  Lastly,  this  material  will 
be  integrated  to  implement  a  quadratic  optimization  MFC  controller. 

2.3.1  Definition  of  MFC  Horizons 

As  a  finite  horizon  control  methodology,  state  space  MFC  requires  the 
introduction  of  three  independent  horizons  in  order  to  function:  the  prediction  horizon, 
the  control  horizon,  and  the  optimization  horizon.  The  prediction  horizon  is  of  length  p 
and  defines  the  number  of  time  steps  into  the  future  for  which  the  plant  state  or  output 
vector  is  predicted  at  any  given  time,  k,  during  the  optimization  process.  The  control 
horizon,  q,  defines  the  number  of  time  steps  over  which  incremental  control  input  moves 
are  available  starting  at  time  k,  after  which  Au(k+l)  =  0  for  /  ^  The  optimization 
horizon,  r,  defines  the  number  of  state  feedback  signals,  v(k+r),  that  are  calculated  by  the 
optimizer.  After  r  steps,  v(^+/)  reaches  the  steady  state  value  v”  that  drives  the  system  to 
its  setpoint,  s”.  The  relative  lengths  of  the  various  horizons  are  important  because  the 
stability  of  the  system  being  controlled  by  the  MFC  process  is  directly  related  to  the 
choice  of  these  three  horizons  [4].  For  a  system  using  output  feedback  and  an  estimator 
for  state  prediction,  stability  of  the  closed  loop  system  operating  under  constraints  is  only 
guaranteed  if:  N  -  2n  <  r  <  min(p-2n,q-2n),  where  n  is  the  order  of  the  plant  and  N  is  the 
fewest  number  of  steps  required  to  drive  the  system  to  its  steady-state  setpoint,  s°. 
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2.3.2  Internally  Stabilizing  Feedback  Controller 


Let  the  discrete  time  state  space  representation  for  a  plant  G{z)  e  be  given  by 

x(A:+l)  =  Ax(k)  +  BAu(k)  n  ^{w 

y{k)  =  Cx{k)  ^  ’ 

with  (A,B)  controllable  and  (C,A)  observable,  and  where  x(k)  e  K"  is  the  system  state 
vector.  In  order  to  implement  the  stabilizing  inner  feedback  controller  described  in 
Section  2.2.2,  we  need  to  find  the  coprime  factorization  of  (2.20).  The  right  coprime 
factorization  may  be  accomplished  through  the  introduction  of  the  state  feedback 
variable,  v(k),  in  the  form  of  equation  (2.10).  F,  is  chosen  such  that  the  poles  of  A  +  BF, 
are  at  the  origin,  which  leads  to  the  recovery  of  the  FIR  relationship  between  the  “quasi¬ 
reference”  input,  v(A:),  and  the  plant  input  and  output.  These  FIR  relationships  allow  for 
the  calculation  of  the  far  future  value  of  v(A:)  that  causes  the  estimated  plant  output  to 
reach  its  setpoint  within  the  prediction  horizon.  The  factors  necessary  to  construct  the 
stabilizing  controller  shown  in  Figure  2.1  are  given  in  equations  (2.15)  and  by 


A  +LC 

-B 

L 

7  * 

0 

(2.21) 


Furthermore,  the  condition  that  Z,Z;'  =  /  is  stipulated  in  order  to  attain  the  closed  loop 
relations  of  equation  (2.14). 
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2.3.3  State  Prediction 


In  order  to  find  the  optimal  set  of  reference  inputs,  v^k+l),  to  attain  and  track  the 
setpoint,  it  is  first  necessary  to  construct  a  methodology  to  determine  the  future  system 
states  based  on  these  inputs.  As  most  physical  systems  have  at  least  some  unmeasurable 
states,  it  is  convenient  to  base  the  prediction  equations  on  an  estimator  of  the  form 

x(k  +  1)  =  Axik)  +  BAu(k)  +  L[Cm  - 
Au(k)  -  F^x(k)  +  Z^v(k) 


with  f  (it)  6  M''.  Substituting  for  Au(k)  in  the  first  equation  in  set  (2.22)  then  leads  to 

xik  +  1)  =  (^  +BF^  +  LC)x(k)  +  BZ^vik)  +  Lyik)  (2.23) 


Using  equation  (2.23),  it  is  now  possible  to  express  the  state  estimate  for  a  particular 
point  /  on  the  prediction  horizon  as 


x(k  +  /  +  1)  =  Fil)x(k  +  /)  +  Gv{k  +  /)  +  H(l)y{k) 

I  =  Q...p- I 


(2.24) 


where  F(0)  =A+BF^  +  LC  and  F(l) ...  F(p-\)  =  A+BF,,G  =  BZ,,  H{0)  =  I  and 
H(\) ...  H(p-\)  =  0.  The  terms  dealing  with  the  output  feedback  are  eliminated  for  all 
times  greater  than  k  because  the  actual  system  output  is  not  available  at  any  future  time. 

Utilizing  equation  (2.24),  it  is  possible  to  construct  a  vector  of  future  state 
estimates  using  the  notation  defined  in  [4]  and  [6],  but  it  is  first  necessary  to  define  the 
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vectors 


v(^)  =  [  v(A:)^  •••  v(k  +  r  -  \Y  Y 
v“(^)  =  [  v(k  +  rY  v{k  +  p  -  1)’"  Y 
Au(k)  =  [  Au(kY  •••  Au(k  +  q  -  \Y  Y 
u(k)  =  [  u{kf  •••  uik  +  q  -  \Y  Y 


(2.25) 


where  p  =  max(p,  ^).  The  vector  of  state  estimate  predictions  can  then  be  written  as 


f(A:  +  1) 

x{k^p) 


Fx{k)  +  Gv{k)  +  G“v“(A:)  +  Hy{k) 


(2.26) 


in  which  F,  G,  G  and  H  are  matrix  functions  of  F(/),  G,  and  H(l).  These  matrix 
functions  can  be  represented  using  the  notation  [4,6] 


Ilfo-) 


j  =  m 


'  F{m)F{m  -  \)...F{n) 
F{m) 

I 

0 


m>  n 
m  -  n 
m  =  n  -  \ 
otherwise 


(2.27) 


which  allows  the  predicted  state  estimate  at  each  time  step  to  be  defined  by 


X(^W  +  1)  =  nF(y)i(^)  +  Si 

7  =  /  /  =  0 


ri^o) 

j=> 


Gv{k  +  i) 


j=‘ 


H(0)y(k)  (2.28) 


Using  equation  (2.28),  it  is  thus  possible  to  formulate  the  state  estimate  prediction 
matrices  F,  G,  G°°,  ^lnd  H: 
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0  0  0 

G  0  0 

F(r  +  l)G  G  ••  0 

G“  =  [  r.l  1  \  P 

nF{j)G  F(r+2)G  •  0 

[j=rt2  J 

n  F(j)  G  n  F(j)  G  •  n  f(j)  g 

J  J  J  J  . 


HiO) 

Fil)HiO) 

H  =  \P 

n  F(j)  H(0) 

J  J  . 


(2.29) 
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Similarly,  it  is  possible  to  construct  a  vector  of  predicted  future  incremental 
control  input  moves,  Au(k) 

Am  =  FJ(k)  +  G^m  +  G;v“  +  n^yik)  (2.30) 

where  F^,  G^,  Gj,  and  are  matrix  functions  of  F^,  F(f),  G,  and  H(l).  To  do  so, 

we  must  first  introduce  a  modified  form  of  the  controller  input  from  equation  (2.22)  that 
is  defined  for  any  time  /  along  the  control  horizon: 

Au{k^l)  =  F^x{k  +  l)  *  Z^v{k^l) 


Inserting  the  state  estimate  prediction  equations  from  (2.24)  into  (2.31)  then  yields  an 
expression  for  the  vector  of  future  incremental  control  inputs 


Au{k  +  /)  =  n  Fo-)m)  +  f.E  1 

y=/-l 


/  =  0 


ri  FU) 

y=/-i 


Z^vik^l)  +  F, 


n  FU) 

j=i-i 


Gv(k  +  i) 

Hi0)yik) 


(2.32) 


from  which  we  may  derive  F^,  G^,  Gj,  and  H^: 


F^  - 


FF{0) 


n  FU) 

j  =  q-2 
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Z 

r 

0 

0 

Ffi 

0 

1 

2 

F 

r 

n  F(/') 

G  F 

n  FU) 

G 

z 

r 

y.r-2 

j^r-2 

1 

2 

Fr 

n  F{j) 

G  F^ 

n  F{j) 

G  •• 

FrG 

i  =  r-\ 

J-r-x 

1 

2 

r 

F. 

n  F{j) 

G  F 

n  FU) 

G  F 

n  FO) 

J  =  q-2 

.V"9-2 

J-q-2 

0 

0 

0 

0 

0 

FrG 

0 

F^F{r  +  \)G 

Ffi 

0 

r  +  1 

r  +  2 

p 

Fr 

n  F{j) 

G  F^ 

n  F<j) 

G  F, 

n  FU) 

J  =  q-2 

(2.33) 
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Lastly,  it  is  possible  to  derive  an  expression  for  the  absolute  control  input  at  any 

/ 

point  /  along  the  control  horizon  based  on  the  fact  that  u(k  + 1)  =  u(k  -  1)  +  '^  Au(k  +j) 

7  =  0 


/  i 

0 

p-i  ( 

\ 

+  0  =  IT  1 

n  Fih) 

m  +  F^ 

n  Fij) 

Gv{k  +  z) 

/;=0  1 

1 

j^h-i 

l.o[ 

[j-h-\ 

Z^v(k  +  h) 

*  F,  n  FU) 

//(0);;(A:)} 

+  u{k  -  1) 

j^h-\ 


(2.34) 


The  vector  of  predicted  absolute  control  inputs  can  thus  be  represented  by 


u(k)  =  F^x(k)  +  G^v(k)  +  Gjv”  +  H^y(k)  +  lu(k-  1)  (2.35) 


where  F^,  G^,  G  and  are  matrix  functions  of  F„  F  (/),  G,  and  H  (/),  and 
/  is  a  column  of  ^  x  ^  identity  matrices; 


Fu  - 


F 


q-\ 

Fr^ 

/  =  0 


n  F{i) 
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FG  +  Z 

r  r 


0 

z 


EF(r) 

(  =  0 


n  F{j) 

7=/-l 


G  +  Zr  EF(r) 
(  =  0 


n  FU) 

j  =  i-\ 


G  +  Zr 


q-\ 


efw 

/  =  0 


n  F{j) 

7  =  <-l 


G  +  Zr  T,F{r) 

i=0 

0 

0 


n  Fij) 

j=i-\ 


G+Zr 


FG  +  Z 

r  r 


q~\ 

LF(r) 

i  =  0 


n  FU) 

j=i-i 


G+Zr 


(2.37) 
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0 


0 


r*\ 

ZF(r) 
/  =  0 


FG.Z 

n  FO) 

y=/-i 


G  +  Zr 


0 

Z 


9-1 

EF(r) 

1=0 


r*\ 

n  F{j) 

7='-l 


G+Zr  EF(r) 
/  =  0 


r*2 

n  F{j) 

j=,-\ 


G  +  Zr 


q-\ 

ZF(r) 

i  =  0 


0 

0 

0 

n  FQ) 

J-'-i 


G+Zr 


H  = 


9-1 

/=0 


0 

FMO) 


n  FU) 

y=/-i 


//(O) 
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2.3.4  System  Constraints 


For  any  physical  system,  there  are  limits  to  the  rate  at  which  a  control  input  may 
be  applied  and  to  the  magnitude  of  the  absolute  control  input.  Additionally,  all  physical 
systems  have  hard  or  soft  constraints  that  may  be  represented  either  by  a  limitation  on  a 
specific  state  or  by  a  limitation  on  a  linear  combination  of  states.  An  example  of  a  hard 
constraint  would  be  the  normal  acceleration  beyond  which  airframe  structural  failure  is 
imminent.  A  soft  constraint  might  be  the  normal  acceleration  that  the  pilot  is  expected  to 
be  able  to  withstand  during  Air  Combat  Maneuvering.  The  rate,  position,  and  state 
estimate  constraints  do  not  necessarily  have  the  same  upper  and  lower  limits  either,  so 
they  are  typically  expressed  as  having  minimum  and  maximum  bounds: 

^  Am(A:  +  i)  <  +  i)  i  =  0...9  -  1 

^  0  ^  +  0  ^  +  0  /  -  0...^  -  1  (2.37) 

4in(^  ^  0  ^  x{k  +  i)  ^  x^^(k  +  0  i  -  1  ...p 

These  bounds  may  be  held  constant  throughout  the  control  and  prediction  horizons,  or 
they  may  be  chosen  to  vary  with  time.  Additionally,  it  is  possible  to  introduce  output 
constraints  on  the  system  by  limiting  linear  combinations  of  the  state  estimate  constraints. 
Furthermore,  if  one  wishes  to  consider  a  particular  control  increment,  control  input,  or 

state  to  be  unconstrained,  the  maximum  or  minimum  limit  may  be  set  to  +0°  or 

respectively. 

The  use  of  the  Matlab  qp  command  for  the  constrained  quadratic  optimization 
process  to  determine  the  series  of  “quasi-reference”  signals,  v(k),  at  each  time  step  makes 
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it  necessary  to  formulate  the  rate,  position,  and  state  constraints  on  the  system  as  a  matrix 
inequality  of  the  form 


Dv{k)  <  Ec{k) 

where  (2.38) 

c{k)  =  [  x(A:)^  y{kY  v'"{kY  u{k-\Y  l{kY  m{kY  n(kY  Y 


and  where  l(k)  E  M^,  m{k)  E  and  n{k)  E  M'’  represent  vectors  of  constraints  on  the 
control  rates,  absolute  control  inputs,  and  state  estimates,  respectively.  Additionally,  X  is 
the  number  of  constraints  on  input  rates  across  the  control  horizon,  p  is  the  number  of 
input  constraints  across  the  control  horizon,  and  v  is  the  number  of  state  estimate 
constraints  across  the  prediction  horizon.  The  elementwise  inequalities  needed  to  obtain 
(2.38)  are 


Lj^Au(k)  <  l(k) 
M^u(k)  <  m{k) 

x(k  +  1) 

■:  <  n(k) 

x(k  +  p) 


(2.39) 


where  in  the  case  without  linear  combinations  of  constraints  L^,  M^,  and  are  defined  by 


(2.40) 
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and  l(k),  m{k),  and  n{k)  are  defined  by: 


m  = 


-U(^) 


,  m{k)  - 


-m  (k) 

mm  ^  ^ 


,  n(k)  = 


min  ^ 


(2.41) 


where  /  e  w  e  and  n  e  Here,  X  and  |a  are  both  assigned  to  be  2q^  because 
there  are  2^  upper  and  lower  input  constraints  at  each  time  step  over  the  control 
horizon,  q.  Similarly,  v  is  set  to  2pK  because  there  are  2k  upper  and  lower  state  estimate 
constraints  at  each  time  step  over  the  prediction  horizon,;?.  By  stacking  the  upper  and 
lower  constraints  as  in  equations  (2.40)  and  (2.41)  and  substituting  the  prediction 
equations  from  (2.26),  (2.30),  and  (2.35)  into  equation  (2.39),  it  is  possible  to  incorporate 
both  the  minimum  and  maximum  bounds  on  the  system  into  the  form  of  equation  (2.38), 
where 


-Ly^H^  -Ly^Gl  0  ''  I  Q  Q 

D  - 

KGu 

,  E  = 

-kg:  -K^  0  ^  0 

N^G 

-N/  -N^H  -N^G^  0  0  0  7 

2.3.5  Quadratic  Programming  Implementation 

In  this  formulation  of  state  space  predictive  control,  the  goal  is  to  find  the  optimal 
series  of  “quasi-reference”  signals,  vik+l),  such  that  the  cost  function 
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(2.43) 


=  E  l|ci(A:  +  /) -5(A:  +  /)||J  +  ^  ||Aw(A:  +  /- 1)||2 

/  =  i  ■'  /=i 

R^>0, 


is  minimized.  R^  and  R^  are  weights  applied  to  the  tracking  and  control  power  terms  of 
the  cost  function  in  order  to  adjust  the  relative  penalties  on  setpoint  tracking  and  control 
usage.  Expanding  this  cost  function  and  substituting  in  the  prediction  equations  (2.26), 
(2.30),  and  (2.35),  one  may  revwite  (2.43)  as  the  equivalent  cost  function 


J(k)  =  v(k/Sv(k)  +  [x(^)^  y(kY  v-(kf  s(kY]Tv(k)  +  K 


(2.44) 


where  S  is  given  by 


5  =  G^C^R^CG  .  GllG, 


(2.45) 


T  is  given  by 


r  =  2 


CF  CH  CG’^  -/ 


p']  J 


^RyCG  +  2 


(?a“  0 


(2.46) 


and  those  terms  not  involving  v(k)  are  grouped  together  in  k: 


K  =  (Fx{k)  +  G“v“  +  Hy(k)yc  ^R^c(Fx(k)  +  G“v”  +  Hy(k))  + 
(FJ(k)  +  g;v“  +  4y(^))X  (FJik)  +  G;v“  +  H^y(k))  + 
5(^)^^Ji^(^)  -2c(Fx(k)  +  Hy(k)  +  G“v“ )] 


(2.47) 
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Additionally,  C  =  diag(C,...,C),  =  diag(i?^,...,i?p,  and  -  diag(i?„,...,i?„). 

Furthermore,  based  on  Lemma  5.1  in  [4],  it  is  possible  to  calculate  the  far  future  value  of 
the  “quasi-reference”  signal,  v”,  from  the  system  of  equations 


In 


CY{A^BFy-'BZ^ 


00  00 

V  =5 


(2.48) 


The  MFC  optimization  may  then  be  realized  as  the  quadratic  programming  problem: 


min  v{kYSv{k)  +  [  y{kY  s{kY  ]Tv(k) 

m 


(2.49) 


subject  to  the  constraints: 


Dv(k)  <  Ec(k) 


(2.50) 


It  should  be  emphasized  at  this  point  that  this  is  not  a  static  quadratic  optimization 
problem.  At  each  time  step,  the  vector  of  state  estimates,  the  plant  output,  the  system 
constraints,  and  the  vectors  of  steady  state  reference  signals  and  setpoints  are  updated  and 
a  new  constrained  quadratic  optimization  is  performed.  These  updates  allow  for  dynamic 
setpoint  and  constraint  changes,  and  the  combination  of  the  on-line  optimization  process 
and  the  stabilizing  inner  feedback  loop  is  the  mechanism  by  which  system  stability  is 
maintained.  The  stable  model  predictive  system  is  implementable  based  on  Figure  2.2. 
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Figure  2.  2  Stable  Model  Predictive  Control  System 
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3.0  System  and  Maneuver  Description  and  Model  Formulation 


3.1  F-18  High  Alpha  Research  Vehicle  ('HARV') 

The  HARV  is  a  NASA-owned  modified  version  of  the  F-18  Hornet  single-seat 
fighter/attack  aircraft  built  by  the  McDonnell  Aircraft  Company,  St.  Louis,  Missouri.  It  is 
powered  by  two  General  Electric  F404-GE-400  afterburning  turbofans,  each  capable  of 
producing  16,000  pounds  static  thrust  in  afterburner  at  sea  level.  In  order  to  augment  yaw 
and  pitch  performance  at  high  angles  of  attack,  three  externally-mounted  thrust  vectoring 
vanes  were  added  to  each  engine  to  deflect  the  exhaust.  The  addition  of  these  paddles 
necessitated  the  removal  of  the  divergent  portion  of  the  engine  nozzles,  which  prevents 
the  HARV  from  achieving  supersonic  flight.  However,  afterburner  operation  is  still 
possible.  A  comparison  of  the  physical  characteristics  of  a  standard  F-1 8  and  the  HARV 
is  included  in  Table  3.1.  A  three- view  drawing  of  the  HARV  may  be  found  in  Figure  3.1. 

The  most  important  functional  improvement  of  the  HARV  over  a  standard  F-1 8  is 
the  addition  of  the  thrust  vectoring  vane  system.  At  low  airspeed,  conventional 
aerodynamic  control  surfaces  are  rendered  ineffective  because  low  dynamic  pressure 
prevents  them  from  generating  the  necessary  moments  to  pitch,  roll,  and  yaw  an  airplane. 
Furthermore,  flight  at  high  angles  of  attack  produces  cross-axis  coupling  of  aerodynamic 
controls,  which  in  turn  complicates  the  development  of  a  flight  control  system.  One 
solution  to  both  these  difficulties  is  to  use  thrust  vectoring.  Thrust  vectoring  can  be  used 
symmetrically  to  produce  either  pitching  moments  or  yawing  moments,  and  it  can  be  used 
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asymmetrically  to  produce  rolling  moments.  Also,  because  the  moments  produced  by 
thrust  vectoring  remain  aligned  with  the  aircraft  axes,  cross-axis  coupling  is  avoided. 
Another  advantage  of  incorporating  thrust  vectoring  into  an  aircraft  is  control 
redundancy.  In  the  event  of  an  actuator  failure  on  one  of  the  aerodynamic  control 
elements,  the  thrust  vectoring  system  is  capable  of  supplementing  the  moments  produced 
by  the  functioning  control  surfaces. 


Table  3.1:  Physical  Characteristics  of  the  Unmodified  and  Modified  F-18 

Parameter 

Unmodified 

Modified 

Weight,  lb 

31,980 

36,099 

Reference  wing  area, 

400 

400 

Reference  m.ax.,  ft 

11.52 

11.52 

Reference  span,  ft 

37.4 

37.4 

Center  of  gravity 

- 

Percent  m.ax. 

21.9 

23.8 

Fuselage  ref.  station 

454.33 

456.88 

Waterline 

105.24 

105.35 

Roll  inertia,  slug-ft^ 

22,040 

22,789 

Pitch  inertia,  slug-ft^ 

124,554 

176,809 

Yaw  inertia,  slug-ft^ 

139,382 

191,744 

Product  of  inertia,  slug-ft^ 

-2,039 

-2,305 

Overall  length,  ft 

56 

56 

Wing  aspect  ratio 

3.5 

3.5 

Stabilator  span,  ft 

21.6 

21.6 

Stabilator  area,  ft^ 

88.26 

86.48 
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Figure  3.1  3-view  drawing  of  the  F-18  HARV 


The  F-18  HARV  is  equipped  with  five  pairs  of  aerodynamic  control  surfaces: 


(1)  Independent  stabilators  capable  of  symmetric  deflection  for  longitudinal  maneuvering 

and  differential  deflection  for  aiding  roll. 

(2)  Single-slotted  ailerons  with  a  ±25  degree  range  of  operation. 
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(3)  Leading-edge  flaps  that  can  deflect  to  a  maximum  of  33  degrees  down.  They 

can  also  operate  differentially  ±3  degrees  for  aiding  roll. 

(4)  Single-slotted  trailing-edge  flaperons  capable  of  a  maximum  45  degree  downward 

deflection  in  the  landing  configuration  and  a  ±8  degree  differential  deflection. 
Additionally,  these  flaperons  are  scheduled  with  Mach  number  and  angle  of  attack 
to  decrease  drag  and  augment  stability. 

(5)  Twin  rudders  used  for  roll  coordination  and  directional  control. 


The  F-18  is  also  equipped  with  large  leading  edge  extensions  (LEX)  that  extend  from  the 
wing  leading  edge  to  just  forward  of  the  cockpit  in  order  to  provide  enhanced 
maneuverability  and  lift  production  at  high  angles  of  attack.  A  complete  listing  of 


position  and  rate  limits  for  the  aerodynamic  control  surfaces  is  presented  in  Table  3.2. 


Table  3.2:  Aerodyriamic  Goiifrol  Surfaces  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

Control  Surface 

Position  Limit,  de& 

Rate  Limit,  deg/sec 

Stabilator 

40 

TEU 

24 

TED 

10.5 

Aileron 

100 

TEU 

25 

TED 

45 

Rudder 

82 

TEL 

30 

TER 

30 

Trailing-edge  flap 

18 

TEU 

8 

TED 

45 

Leading-edge  flap 

15 

LEU 

3 

LED 

33 
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The  non-aerodynamic  controls  of  the  F- 1 8  H  ARV  are  throttle  position  and  thrust 
vectoring  angle.  Information  on  these  control  elements  is  provided  in  Table  3.3. 


Table  3.3:  Propulsive  Control  Elements 

Thrust  Vectoring  Vanes 

Dimension,  in. 

Upper 

20  by  20 

Lower 

20  by  15 

Area,  in^ 

Upper 

358.76 

Lower 

263.64 

Position  Limit,  deg 

Upper 

-10 

Lower 

25 

Rate  Limit,  deg/sec 

80 

Throttle 

Position  Limit,  deg 

Upper 

127 

Lower 

54 

Rate  Limit,  deg/sec 

30 

12  F-18HARV  Model 

Summing  forces  and  moments  acting  on  an  aircraft,  the  non-linear  equations  of 
motion  can  in  general  be  represented  by 


xiO  =/(x(0,m(0) 


(3.1) 


where  x{t)  is  the  vector  of  aircraft  states  and  u(f)  is  the  vector  of  aircraft  control  inputs. 
Furthermore,  at  any  arbitrary  equilibrium  point  (x^,  wj,  the  derivative  of  the  state  vector  is 
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equal  to  the  zero  vector: 


=  0  (3.2) 

Assuming  small  perturbations  to  the  aircraft  states  and  taking  equation  (3.2)  into  account, 
we  can  use  a  Taylor  series  to  expand  equation  (3.1)  and,  neglecting  higher  order  terms  in 
the  expansion,  linearize  the  aircraft  equations  of  motion  about  an  equilibrium  point.  It  is 
then  possible  to  represent  small  motions  about  this  equilibrium  point  with  the  state  space 
system 

x(t)  -  Ax{t)  +  Bu(t) 
y(t)  =  Cx(f)  +  Du{t) 

where  x{t)  and  u{t)  now  represent  perturbations  in  the  states  and  inputs  about  their  trim 
values. 

The  F-1 8  HARV  state  space  models  used  in  this  research  effort  are  taken  from 
[7],  but  were  originally  provided  by  NASA  Langley.  These  models  represent  both 
longitudinal  and  lateral-directional  aircraft  behavior  about  three  trim  conditions  described 
in  Section  3.4.  In  these  state  space  systems,  ^  e  M 5  e  M the  vector  of  aircraft 
states  is 

X  -  ^Vj.  a  P  p  q  r  ^  0  h  Y  (3.4) 

and  the  vector  of  aircraft  inputs  is 

“  "  [  ^TVL  ^TVR  ^RL  ^RR  ^AL  ^AR  ^SL  ^SR  ^LEL  ^LER  ^TEL  ^TER  ]  (3*5) 
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Table  3.4  contains  definitions  of  these  quantities  and  their  associated  units. 


Table  3.4:  Definitions  of  F-18  HARV  States  and  Control  Inputs 

States: 

Vr 

Perturbation  in  true  airspeed 

ft/sec 

a 

Perturbation  in  angle  of  attack 

rad 

p 

Perturbation  in  sideslip  angle 

rad 

p 

Perturbation  in  roll  rate 

rad/sec 

q 

Perturbation  in  pitch  rate 

rad/sec 

r 

Perturbation  in  yaw  rate 

rad/sec 

4> 

Perturbation  in  roll  angle 

rad 

e 

Perturbation  in  pitch  angle 

rad 

Perturbaion  in  yaw  angle 

rad 

h 

Perturbation  in  altitude 

ft 

1  Controls:  | 

^TVL 

Perturbation  in  left  thrust  vectoring  vane  deflection 

deg 

^TVR 

Perturbation  in  right  thrust  vectoring  vane  deflection 

deg 

^Rl. 

Perturbation  in  left  rudder  deflection 

deg 

Perturbation  in  right  rudder  deflection 

deg 

Perturbation  in  left  aileron  deflection 

deg 

^AR 

Perturbation  in  right  aileron  deflection 

deg 

^SL 

Perturbation  in  left  stabilator  deflection 

deg 

^SR 

Perturbation  in  right  stabilator  deflection 

deg 

^IML 

Perturbation  in  left  leading  edge  flap  deflection 

deg 

^LER 

Perturbation  in  right  leading  edge  flap  deflection 

deg 

^TEL 

Perturbation  in  left  trailing  edge  flap  deflection 

deg 

^TER 

Perturbation  in  right  trailing  edge  flap  deflection 

deg 

Perturbation  in  throttle  position 

deg 
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For  a  conventional  aircraft  in  which  each  control  surface  pair  functions  as  a  single 
entity,  linearization  of  the  aircraft  equations  of  motion  about  a  trim  condition  permits 
the  assumption  that  the  longitudinal  and  lateral-directional  dynamics  are  decoupled. 
However,  in  the  case  of  an  aircraft  such  as  the  F-1 8  HARV,  each  control  surface  element 
is  capable  of  operation  independent  of  its  counterpart  and  when  deflected  excites  both 
longitudinal  and  lateral  dynamics.  In  order  to  treat  the  longitudinal  and  lateral-directional 
dynamics  of  the  HARV  independently  after  linearization,  it  is  necessary  to  decompose  the 
motion  of  each  pair  of  control  elements  into  its  S)mimetric  and  differential  components. 
The  symmetric  component  affects  longitudinal  dynamics,  whereas  the  differential 
component  produces  lateral-directional  motion.  Define  these  components  as 


where  X  represents  an  arbitrary  control  element,  R  and  L  specify  the  right  or  left  control 
surface,  D  indicates  differential  deflection,  and  S  indicates  symmetric  deflection.  Using 
these  definitions,  it  is  then  possible  to  formulate  the  decoupled  state  space  system 

^LONO  ^LONG  ^  ^LONG  ^LONG  ^  ^LONG 

—  + 

^LAT  ^  -^UT  ^LAT  ^  ^LAT  ^LAT 


where  e  ^  e  the  longitudinal  state  vector  and 


^LONG  [  CC  0  /?  j  ^LONG  [  ^TVS  ^SS  ^LES  ^TES  ^T  ]^  (3*8) 

and  the  lateral-directional  state  vector  and  input  vector  are 

^LAT  “  [  P  P  ^  ]  ^lAT  ~  [  ^rVD  ^AD  ^SD  ^LED  ^TED  ^RD  ^RS  (3*9) 

Only  longitudinal  dynamics  are  considered  in  this  thesis,  so  extracting  the  longitudinal 
portion  of  the  plant  and  input  matrices  of  the  system  of  equation  (3.7)  and  ignoring  the 
redundant  state,  h,  yields 

y(‘)  = 

where  .4  e  K**  5  e  M'*  ’‘®,  C  e  the  state  and  input  vectors  are 

^LONG  “  [  ^  9  0  ]  ^LONG  ~  [  ^TVS  ^AS  ^SS  ^LES  ^TES  ]  (3-11) 

and  the  output  vector  is  chosen  to  be 

y  =[V^  y  e  (3.12) 

where  (y  =  0  -  a)  is  the  perturbation  in  the  flight  path  angle.  The  output  vector  shown  in 
equation  (3.12)  is  chosen  in  order  to  facilitate  the  establishment  of  setpoints  required  to 
accomplish  the  longitudinal  maneuvers  employed  in  this  thesis.  This  choice  of  outputs 
also  ensures  that  ( C,  ^ ) ,  the  discrete  time  output  and  plant  matrices  of  the  system 
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modified  to  accept  control  increments,  is  observable  for  the  flight  conditions  used  in  this 
thesis.  It  should  be  noted  that  these  perturbed  quantities  are  not  necessarily  directly 
measurable  in  reality;  however,  in  this  thesis  they  are  assumed  to  be  measurable  and 
available  for  feedback. 

3.3  Input-Output  and  State  Scaling 

When  trying  to  establish  the  relative  importance  of  output  deviations  from  a 
setpoint,  it  is  beneficial  to  employ  units  or  non-dimensionalized  quantities  that  permit  an 
accurate  comparison  of  these  deviations.  For  instance,  a  1  ft/sec  deviation  in  the  forward 
velocity  of  an  aircraft  is  certainly  not  comparable  to  a  1  radian  deviation  in  pitch  angle. 
Additionally,  given  the  substantial  differences  in  the  absolute  position  limits  of  the 
control  surfaces,  scaling  the  control  deflections  based  on  these  limits  or  on  other  more 
convenient  quantities  allows  better  visualization  of  relative  control  usage.  After  the 
appropriate  scale  factors  are  chosen  with  regard  to  a  particular  trimmed  flight  condition, 
the  input-output  and  state  scalings  may  be  accomplished  through  elementary  coordinate 
transformation  procedures.  Consider  the  state  space  realization 

(t)  + 

'  (3-13) 

with  original  state,  input,  and  output  quantities  (x^,  m„,  Their  scaled  coimterparts 
(x,  y,  u)  are  related  to  the  original  quantities  through  the  relations 
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(3.14) 


=  Tx 
u  =  Su 

O 


where  T,  S,  and  R  are  transformation  matrices  of  appropriate  dimensions.  Substituting  the 
transformation  equations  from  (3.14)  into  equation  (3.13)  yields  the  following  scaled 
system 


x(f)  =  Ax(t)  +  Bu{t) 
y{f)  =  Cx{t) 


(3.15) 


where 


A  =  T-Uj  B  =  T^^B^S  C  =  R-^CJ  (3.16) 


In  this  thesis,  the  state  transformation  from  radians  to  degrees  was  the  first  step  in 
the  scaling  process.  Next,  with  the  exception  of  the  leading  edge  and  trailing  edge  flaps, 
the  inputs  were  scaled  based  on  twice  the  absolute  difference  between  the  trim  deflection 
and  the  nearest  absolute  position  limit.  The  factor  of  two  comes  from  the  definition  of 
the  symmetric  control  deflection  in  equation  (3.6).  The  leading  edge  and  trailing  edge 
flaps  are  special  cases  because,  although  the  nearest  absolute  position  limits  are  upward 
(negative)  deflections,  the  movement  of  these  control  surfaces  during  the  commanded 
maneuvers  is  primarily  in  the  downward  (positive)  direction.  Therefore,  scaling  factors 
for  these  two  control  surfaces  were  chosen  to  enhance  the  clarity  of  the  control  signal 
graphs  and  then  multiplied  by  a  factor  of  two  for  consistency  with  the  other  aerodynamic 
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control  surfaces.  Finally,  the  outputs  were  scaled  to  depict  the  deviations  of  the  flight 
path  angle  and  the  pitch  angle  from  the  trim  condition  in  degrees,  and  the  scaling  of  the 
deviation  in  forward  velocity  from  trim  was  determined  based  on  an  estimation  of  its 
relative  importance  with  respect  to  a  deviation  in  the  flight  path  and  pitch  angles.  A 
complete  listing  of  all  the  scaling  factors  for  the  flight  condition  used  in  this  thesis  is 
contained  in  Table  3.5. 

One  consequence  of  input-output  scaling  that  must  be  taken  into  account  is  its 
effect  on  the  performeince  index  used  in  the  quadratic  optimization  procedure  delineated 
in  Section  2.3.6.  Recall  the  performance  index,  J{k),  is  defined  by 

J{k)  =  £  |cf(^  +  /)-5(^  +  /)||^  +  £  ||Am(^W-1)||^  (3.17) 

/=i  ^  i=\ 

Now,  note  that  each  term  in  the  summation  of  the  control  increment  inputs  over  the 
control  horizon  is  of  the  form 

Lu{k  +  tf  R^Lu(k  +  l)  (3.18) 

where  Aw  refers  to  the  scaled  control  increment  inputs  and  is  a  weight  on  control 
usage.  Substituting  the  control  scaling,  Aw^  =  5Aw,  in  which  Aw^  refers  to  the  unsealed 
control  input  increments,  into  equation  (3.18)  then  yields  the  modified  cost  term 

^u^{k  +  +  1)  (3.19) 
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which  clearly  indicates  a  new  weighting  factor,  R  =  \  is  being  applied  to  the 

unsealed  control  increment  inputs.  Thus,  the  scaling  factors  distinctly  influence  the 
distribution  of  control  power  to  those  control  surfaces  with  the  largest  scaling  factors.  To 
rectify  this  situation,  the  desired  control  power  weighting  should  be  established  in  R,  and 
should  then  be  calculated  by 

R^=S^RS  (3.20) 

Unless  otherwise  indicated,  R  will  always  be  chosen  as  a  constant,  c,  times  the  identity 
matrix  in  this  thesis,  which  further  simplifies  equation  (3.20)  to 

R^  =  cS'^S  (3.21) 

The  scaled  and  unsealed  models  may  be  found  in  Appendices  A1  and  A2. 


Table  3.5:  Scaling  Factors 

Variables 

Fit.  Cond.  1 

Vr 

1/8  ft/sec 

Y 

1/1  deg 

e 

1/1  deg 

6tvs 

1/20  deg 

^AS 

1/50  deg 

^ss 

1/34  deg 

^LES 

1/30  deg 

^TES 

1/60  deg 

bj 

1/27  deg 
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3.4  Flight  Condition 


The  flight  condition  utilized  in  this  research  effort  represents  level  flight  at  an 
altitude  of  15,000  feet  at  Mach  0.24.  It  is  expected  to  be  a  challenging  case  because  the 
low  airspeed  and  the  high  angle  of  attack  diminish  the  effectiveness  of  the  aerodyneimic 
control  surfaces  and  also  because  the  F-1 8  HARV  has  an  unstable  phugoid  mode  at  this 
trim  point.  As  a  consequence  of  the  decreased  aerodynamic  control  authority,  thrust 
vectoring  should  play  an  important  role  in  performing  the  selected  maneuvers.  Table  3.6 
contains  the  initial  aircraft  states  and  control  deflections  for  this  trim  point,  and  Table  3.7 
lists  its  open  loop  poles,  damping  ratio,  and  natural  frequency. 


Table3.6:  iTnm  l 

(lighif  Conditions  ' 

Parameter 

Altitude  (ft) 

15000 

Mach  Number 

0.24 

Vj  (ft/sec) 

238.7 

a (deg) 

25 

q  (deg/sec) 

0 

e  (deg) 

25 

Y  (deg) 

0 

brvs  ~  ^TVL  ^TVR 

0  +  0 

^AS  ~  ^AL  ^AR 

0  +  0 

^SS  ~  ^SL  ^SR 

(-6.40)  +  (-6.40) 

^lES  ~  ^LEL  ^  ^LER 

(0.33)  +  (0.33) 

^TES  ~  ^TEL  ^TER 

(0.80)  +  (0.80) 

100.5 
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Table  3.7:  Open  Loopholes 

Parameter 

Phugoid  poles 

0.0188  ±0.1280j 

Phugoid  damping 

N/A 

Phugoid  nat.  freq. 

0.13 

Short  period  poles 

-0.2481  ±0.3585j 

Short  period  damping 

0.57 

Short  period  nat.  freq. 

0.44 

3.5  Precision  Longitudinal  Maneuvers 

The  three  longitudinal  precision  control  modes  of  interest  to  fighter  aircraft 
employed  in  this  thesis  are  vertical  translation,  pitch  pointing,  and  direct  lift,  of  which 
descriptions  are  taken  from  [7].  During  a  vertical  translation  maneuver,  the  pitch  angle  is 
held  constant  while  the  flight  path  angle  is  varied.  Pitch  pointing  is  characterized  by  a 
change  in  pitch  angle  with  no  variation  in  flight  path  angle.  Lastly,  direct  lift  is  a  mode  in 
which  the  flight  path  angle  and  the  pitch  angle  change  so  that  angle  of  attack  is  held 
constant.  The  setpoints  based  on  the  outputs  specified  in  equation  (3.12)  may  be  found  in 
Table  3.8. 


Table  3.8:  Setpoint  Specifications 

Vertical  Translation 

=  [  0  1  0 

Pitch  Pointing 

11 

O 

O 

H 

Direct  Lift 

=  [  0  1  1 
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4.0  Simulation  Results 


The  following  sections  present  the  results  of  simulations  performed  in  order  to 
determine  the  effectiveness  of  a  state  space  MFC  controller  at  handling  actuator  failures. 
All  six  available  control  inputs  are  used  in  every  simulation  because  this  scenario  offers 
the  greatest  amount  of  control  redimdancy,  which  is  especially  important  in  the  case  of 
out-of-trim  failure  conditions.  To  simulate  these  failures,  multiple  sets  of  constraints  are 
formulated,  and  at  the  time  of  the  simulated  failure,  the  optimization  function  is  updated 
to  operate  using  the  constraint  set  that  defines  the  failure  condition.  Additionally,  if 
necessary  to  improve  performance  in  the  post-failure  operating  environment,  the 
weighting  matrices  Ry  and  may  be  updated  in  the  optimization  through  the  introduction 
of  new  S  and  T  terms  as  shown  in  equations  (2.45)  and  (2.46). 

Failures  are  typically,  but  not  always,  specified  by  setting  the  rate  limit  on  a 
particular  control  element  to  a  very  small,  but  non-zero  value.  In  this  fashion,  the  failures 
are  made  to  occur  at  natural  points  along  the  travel  of  the  control  surfaces,  and  artificially 
extreme  or  conservative  failure  cases  are  avoided.  In  general,  single  element  failures  are 
specified  at  the  maximum  travel  of  the  specific  element  during  a  longitudinal  maneuver. 
The  time  and  magnitude  of  multiple  element  failures  are  determined  on  a  case  by  case 
basis. 

Because  of  the  results  in  [4],  constraints  are  applied  only  at  the  current  time,  k. 
Application  of  constraints  across  the  entire  prediction  and  control  horizons  often  leads  to 
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a  more  conservative  approach  to  the  tracking  problem,  which  is  undesirable  with  regards 
to  a  fighter-type  aircraft.  Additionally,  when  quadratic  programming  is  used  as  the 
method  for  finding  the  series  of  quasi-reference  signal  inputs,  infeasibility  of  the 
optimization  results  if  the  system  is  over-constrained.  For  the  purposes  of  this  thesis, 
feasibility  is  a  precondition  imposed  on  all  the  simulations  presented  here.  For 
information  regarding  the  handling  of  infeasibility,  see  chapter  6  in  [4]. 

MFC  with  quadratic  programming  as  the  optimization  method  is  veiy  susceptible 
to  problems  related  to  numerical  conditioning  and  constraint  application.  Overly 
stringent  constraints  will  prevent  the  optimization  from  finding  a  feasible  solution  and  the 
system  will  consequently  become  unstable.  Therefore,  when  failures  are  specified,  care 
should  be  taken  to  not  set  the  upper  and  lower  rate  limits  to  zero,  for  example.  Instead, 
rate  limits  of  1 0  to  1 0  *  are  sufficient  to  make  the  control  surface  ineffective  and 
generally  do  not  lead  to  numerical  instability.  Furthermore,  the  optimization  weightings 
should  not  be  chosen  such  that  the  quadratic  programming  matrices  are  numerically  ill 
conditioned.  The  results  of  the  simulations  will  be  unreliable  and  instability  often  occurs. 

Many  simulations  were  performed  for  this  thesis  in  which  the  attempt  was  made 
to  decrease  the  sample  time  of  the  system  and  thus  improve  overall  system  tracking 
performance.  In  all  cases,  decreasing  the  sample  rate  below  0.5  seconds  did  not  improve 
performance  and  very  often  resulted  in  worse  performance.  The  reason  is  primarily  that 
the  pole  placement  algorithm  becomes  less  and  less  efficient  as  the  sample  time  is 
decreased,  causing  the  poles  to  be  further  from  the  origin  than  is  desirable.  Hence,  FIR 
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characteristics  are  diminished,  and  tracking  performance  suffers.  Another  result  of  the 
lower  sample  time  is  to  make  the  system  very  susceptible  to  changes  in  the  optimization 
weightings.  Small  changes  in  the  weights  often  lead  to  large  overshoots  and/or  non¬ 
minimum  phase  behavior. 

All  relevant  Matlab  functions  are  contained  in  Appendices  B1  through  B8.  A 
diagram  of  the  simulink  model  used  to  perform  the  simulations  is  contained  in 
Appendix  C. 

4.1  Vertical  Translation 

As  mentioned  earlier,  vertical  translation  represents  a  change  in  the  flight  path 
angle  of  the  aircraft  while  holding  the  forward  velocity  and  the  pitch  angle  constant.  The 
added  lift  production  is  primarily  from  the  ailerons  and  the  trailing  edge  flaps,  with  all  the 
constrained  cases  favoring  the  ailerons  because  of  their  more  responsive  rate  limit.  The 
results  of  [4,7]  and  the  simulations  in  this  thesis  indicate  that  thrust  vectoring  is  used 
primarily  to  maintain  the  proper  pitch  angle  and  minimize  deviations  of  the  forward 
velocity  from  the  trim  condition.  In  all  vertical  translation  cases,  Ry  =  diag(50,20,10''), 

=  (lO'^)S^S,  and  the  sample  rate  is  0.5  seconds.  The  small  control  weightings  increase 
the  overall  speed  of  the  system  response  and  often  cause  the  control  surfaces  to  reach 
their  rate  and/or  position  limits  during  the  aggressive  tracking  behavior  that  results.  The 
commanded  input  is  a  unit  pulse  with  a  duration  of  five  seconds. 

Figures  4. 1  through  4.4  represent  a  case  in  which  the  optimization  is  performed 
without  imposing  any  saturation  or  rate  limits  on  the  control  surfaces.  The  aircraft 
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Time  (sec) 

Figure  4.1  VT  Output  Response  (Unconstrained) 


Time  (sec) 

Figure  4.2  VT  Control  Deflections  (Unconstrained) 


Time  (sec) 

Figure  4.3  VT  Control  Deflections  (Unconstrained) 


Time  (sec) 

Figure  4.4  VT  Control  Deflections  (Unconstrained) 


reaches  its  setpoint  within  0.5  seconds,  but  not  without  violating  rate  limits  on  both  the 
throttle  and  the  trailing  edge  flaps  and  the  upper  saturation  limit  on  the  trailing  edge  flaps. 
The  application  of  the  nominal  constraint  set,  as  shown  in  Figures  4.5  through  4.8, 
produces  a  large  increase  in  overall  control  usage.  Because  the  flaps  are  now  heavily 
constrained  by  their  rate  limit,  the  ailerons  are  forced  to  their  maximum  downward 
position  limit  in  order  to  produce  the  lift  required  to  accomplish  the  maneuver. 
Furthermore,  during  the  initial  maneuver  to  reach  the  setpoint,  rate  limits  are  encountered 
in  both  the  trailing  edge  flaps  and  the  throttle,  and  the  thrust  vectoring  vane  attains  its 
upward  position  limit  in  order  to  counteract  the  effects  of  the  negative  pitching  moment 
produced  by  the  downward  deflections  of  the  ailerons  and  flaps.  Similarly,  when  the 
aircraft  is  commanded  to  return  to  its  original  equilibrium  condition  at  time  t  =  5  seconds, 
the  ailerons  are  driven  to  their  maximum  upward  position  limits  in  order  to  dissipate  lift, 
while  the  thrust  vectoring  vanes  attain  their  maximum  downward  limit  to  counteract  the 
positive  pitching  moment  generated  by  the  flaps  and  ailerons  and  to  prevent  large 
deviations  in  forward  velocity. 

The  next  three  cases  involve  failures  of  single  control  surfaces.  In  Figures  4.9 
through  4.12,  the  stabilator  is  frozen  at  time  t  =  0.5  seconds.  The  flaps  respond  by 
returning  to  a  position  near  their  original  trim  point  to  diminish  the  effects  of  the  extra  lift 
caused  by  the  out-of-trim  failure  of  the  ailerons.  Additionally,  the  throttle  is  retarded  to 
help  prevent  the  flight  path  angle  from  exceeding  its  setpoint.  The  stabilator  and  the 
thrust  vectoring  vane  are  pushed  almost  to  their  upper  and  lower  limits,  respectively, 
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Scaled 


Figure  4.8  VT  Control  Deflections  (Nominal  Constraints) 
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Scaled  Inputs 


Time  (sec) 

Figure  4.9  VT  Output  Response  (Failure:  AS  at  t=0.5) 


Figure  4.10  Control  Deflections  (Failure:  AS  at  t=0.5) 
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Figure  4.12  VT  Control  Deflections  (Failure:  AS  at  t=0.5) 


subsequent  to  the  aileron  failure.  Most  likely,  the  stabilator  is  attempting  to  counteract 
the  effects  of  the  negative  pitching  moment  created  hy  the  downward  deflection  of  the 
ailerons  while  the  thrust  vectoring  vanes  attempt  to  minimize  deviations  of  the  forward 
velocity  from  trim.  This  failure  has  little  effect  on  the  initial  setpoint  tracking  other  than  a 
slight  overshoot  that  does  not  occur  in  the  nominal  case.  However,  the  return  to  the 
original  trim  condition  takes  almost  1.5  seconds  longer  due  to  the  fact  that  the  flaps  must 
work  against  their  rate  limit  to  dissipate  the  excess  lift. 

Figures  4.13  through  4.16  present  a  failure  scenario  in  which  the  stabilator  fails  at 
time  t  =  1 .0  second.  The  relatively  small  deflection  at  which  the  stabilator  freezes  is 
easily  counteracted  by  the  thrust  vectoring  vane  system,  and  the  only  noticeable  effects  on 
the  output  response  are  a  slightly  longer  time  required  to  return  to  the  original  trim 
condition  and  a  greater  deviation  from  the  trim  forward  velocity  than  in  the  nominal  case. 

The  last  single  failure  case  presented  here  is  a  failure  of  the  trailing  edge  flaps  at 
time  t  =  1 .0  second:  Figures  4. 1 7  through  4.20.  Of  those  simulations  included  in  this 
thesis,  this  scenario  presents  the  most  difficult  case  for  the  MFC  controller  in  terms  of 
single  control  element  malfunctions.  In  response  to  the  update  of  the  controller  with  the 
new  constraint  set,  the  ailerons  are  immediately  deflected  upward  in  order  to  shed  lift  and 
prevent  the  aircraft  from  overshooting  its  setpoint.  Simultaneously,  the  stabilator  deflects 
upward  to  counteract  the  negative  pitching  moment  produced  by  the  flap  deflection,  and 
the  thrust  vectoring  vane  closely  approaches  its  lower  limit  to  minimize  any  airspeed 
deviation.  From  the  output  response,  we  see  that  there  is  a  slight  overshoot,  but  the 
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Figure  4.16  VT  Control  Deflections  (Failure:  SS  at  t=1.0) 
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Figure  4.17  VT  Output  Response  (Failure:  TES  at  t=1.0) 
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Figure  4.18  VT  Control  Deflections  (Failure:  TES  at  t=1.0) 
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Figure  4.19  VT  Control  Deflections  (Failure:  TES  at  t=1.0) 
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Figure  4.20  VT  Control  Deflections  (Failure:  TES  at  t=1.0) 


airspeed  deviation  from  trim  is  not  significantly  worse  than  that  of  the  nominal  case.  The 
trouble  arises  when  it  is  time  to  return  to  the  original  trim  condition.  The  ailerons, 
already  near  their  maximum  upward  deflection,  are  pushed  completely  to  saturation  in 
order  to  dissipate  enough  lift  to  return  the  flight  path  angle  to  horizontal.  The  throttle  is 
also  retarded  beyond  its  position  in  the  nominally  constrained  case  because  less  power  is 
required  to  maintain  level  flight  with  the  flaps  deflected. 

Now  that  we  have  shown  the  MFC  controller  is  capable  of  handling  single  control 
element  malfunctions  with  minimal  performance  loss,  it  is  time  to  explore  multiple 
failure  scenarios.  The  first  of  these  is  shown  in  Figures  4.21  through  4.24.  In  this 
instance,  the  stabilator  is  assumed  to  malfunction  at  time  t  =  0  seconds  and  is 
immediately  followed  by  an  out-of-trim  failure  of  the  trailing  edge  flaps  at  t  =  0.5 
seconds.  An  interesting  result  of  a  failure  of  the  stabilator  at  the  trim  condition  prior  to  a 
vertical  translation  maneuver  is  an  implied  constraint  on  the  movement  of  the  thrust 
vectoring  vanes.  Considering  the  large  weight  placed  on  deviations  of  the  pitch  angle 
from  trim  in  these  simulations,  the  controller  carmot  move  the  thrust  vectoring  vanes  to 
any  great  degree  without  significantly  altering  the  pitch  angle.  In  a  similar  simulation  not 
included  here,  a  failure  of  the  thrust  vectoring  system  prior  to  the  maneuver  produces  an 
implied  constraint  on  stabilator  movement  as  well.  From  the  output  response,  we  see  that 
the  primary  effects  of  this  failure  are  similar  to  those  of  the  case  where  only  the  trailing 
edge  flaps  fail.  This  scenario  does,  however,  have  a  longer  rise  time  and  a  slightly  larger 
deviation  in  forward  velocity  from  trim  during  the  first  one  second  of  the  simulation. 
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Scaled  Inputs 
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Figure  4.21  VT  Output  Resp.  (Failures:  SS  at  t=0;  TES  at  t=0.5) 
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Figure  4.22  VT  Control  Defl.  (Failures:  SS  at  t=0;  TES  at  t=0.5) 
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Figure  4.23  VT  Control  Defl.  (Failures:  SS  at  t=0;  TES  at  t=0.5) 
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Figure  4.24  VT  Control  Defl.  (Failures:  SS  at  t=0;  TES  at  t=0.5) 
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The  next  dual  failure  condition,  shown  in  Figures  4.25  through  4.28,  presents  the 
interesting  case  of  a  complete  loss  of  primary  pitch  control  at  the  trim  position.  At  time 
t  =  0  seconds,  the  thrust  vectoring  vanes  and  the  stabilator  are  frozen  and  the  aircraft  is 
commanded  to  perform  vertical  translation.  Because  vertical  translation  does  not  involve 
a  pitch  change,  the  aircraft  is  still  capable  of  completing  the  maneuver  and  returning  to  its 
original  equilibrium  point.  From  Figures  4.27  and  4.28,  we  observe  that,  compared  with 
the  nominal  case,  the  maximum  aileron  deflection  is  reduced  whereas  the  maximum  flap 
deflection  is  increased.  Presumably  this  change  in  lift  production  distribution  occurs  in 
order  to  diminish  the  negative  pitching  moment  produced  by  downward  deflections  of  the 
flaps  and  ailerons.  The  upward  deflection  of  the  leading  edge  flaps  appears  to  counteract 
the  negative  pitching  moment  that  is  produced  by  the  trailing  edge  flaps  and  the  ailerons. 

The  scenario  presented  in  Figures  4.29  through  4.32  is  similar  to  the  previous  case 
except  that  in  addition  to  a  complete  primary  pitch  control  failure,  the  trailing  edge  flaps 
are  frozen  at  time  t  =  0.5  seconds.  From  the  output  response  in  Figure  4.29,  it  is  apparent 
that  this  series  of  failures  presents  a  serious  hindrance  to  nominal  performance.  More 
importantly,  however,  the  aircraft  remains  stable  and  the  controller  is  eventually  able  to 
return  the  aircraft  to  level  flight.  Control  usage  is  similar  to  that  of  the  preceding  failure 
scenario,  but  a  lower  throttle  setting  and  a  saturated  upward  aileron  deflection  are 
necessary  to  return  the  flight  path  angle  to  zero. 

The  last  three  failure  case  is  displayed  in  Figures  4.33  through  4.36:  the  stabilator 
fails  at  time  t  =  0;  the  leading  edge  flaps  fail  at  time  t  =  0.5;  the  ailerons  fail  at  t  =  1.0. 
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Time  (sec) 

Figure  4.25  VT  Output  Resp.  (Failures:  SS  and  TVS  at  t=0) 


Time  (sec) 

Figure  4.26  VT  Control  Defl.  (Failures:  SS  and  TVS  at  t=0) 


Time  (sec) 

Figure  4.27  VT  Control  Defl.  (Failures:  SS  and  TVS  at  t=0) 


Scaled  Inputs 


Figure  4.32  Controls  (Failures:  SS  and  TVS  at  t=0;  TES  at  t=0.5) 
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Time  (sec) 

Figure  4.33  Outputs  (Failures:  SS  at  t=0;  LES  at  t=0.5;  AS  at  t=1.0) 
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Time  (sec) 


Figure  4.35  Controls  (Failures:  SS  at  t=0;  LES  at  t=0.5;  AS  at  t=1.0) 
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Time  (sec) 

Figure  4.36  Controls  (Failures:  SS  at  t=0;  LES  at  t=0.5;  AS  at  t=1.0) 
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Studying  the  control  inputs  in  Figures  4.34  through  4.36,  we  once  again  see  the  implied 
limitation  on  the  deflections  of  the  thrust  vectoring  vanes  because  the  stabilator  is  frozen 
in  its  trim  position  and  also  because  of  the  high  weighting  term  on  deviations  in  pitch 
angle.  However,  because  of  the  negative  pitching  moment  produced  by  the  failure  of  the 
ailerons  in  an  out-of-trim  position,  the  thrust  vectoring  vanes  are  canted  upward 
throughout  the  simulation.  This  movement  limitation  leads  to  significantly  increased 
deviations  in  forward  velocity  from  the  nominally  constrained  case  ^lnd  contributes  to  the 
overshoot  shown  in  the  output  response  in  Figure  4.33.  Additionally,  the  failed  ailerons 
are  responsible  for  the  slow  return  to  the  original  equilibrium  condition. 

Except  for  the  single  control  element  failures,  most  of  the  failure  scenarios 
presented  in  this  thesis  are  highly  unlikely  unless  the  aircraft  suffers  some  form  of 
external  damage.  However,  the  benefit  from  exploring  these  unlikely  cases  is  derived 
from  the  fact  that  a  controller  capable  of  handling  such  severe  failures  should  be  readily 
capable  of  handling  minor  ones.  Still,  it  is  worthwhile  to  explore  some  instances  in 
which  the  MFC  controller  is  either  unable  to  return  the  aircraft  to  its  original  trim  point  or 
is  unable  to  prevent  the  aircraft  from  departing  from  controlled  flight.  The  first  of  these 
cases  is  shown  in  Figures  4.37  through  4.40  and  represents  a  failure  of  the  ailerons  at 
t  =  0.5  seconds  and  a  malfunction  of  the  trailing  edge  flaps  at  t  =  1.0  second.  Noting  that 
these  two  pairs  of  control  surfaces  are  those  primarily  responsible  for  generating  the  extra 
lift  required  for  the  vertical  translation  maneuver,  we  should  expect  that  both  tracking  the 
setpoint  and  returning  to  equilibrium  should  be  nearly  impossible,  and  this  is  the  case. 
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Time  (sec) 

Figure  4.37  VT  Output  Resp.  (Failures:  AS  at  t=0.5;  TES  at  t=1.0) 


Time  (sec) 


Figure  4.38  VT  Control  Defl.  (Failures:  AS  at  t=0.5;  TES  at  t=1.0) 


Scaled  Inputs 


Time  (sec) 

Figure  4.39  VT  Control  Defl.  (Failures:  AS  at  t=0.5;  TES  at  t=1.0) 


Time  (sec) 


Figure  4.40  VT  Control  Defl.  (Failures:  AS  at  t=0.5;  RES  at  t=1.0) 


Subsequent  to  the  trailing  edge  flap  failure  at  time  t  ==  1 .0  second,  the  stabilator,  the  thrust 
vectoring  vanes,  and  the  leading  edge  flaps  all  reach  either  their  positive  or  negative 
saturation  limits.  As  seen  in  the  output  response  in  Figure  4.37,  even  these  extreme 
deflections  are  ineffective  when  it  comes  to  tracking  the  setpoint  or  returning  to 
equilibrium.  The  case  shown  in  Figures  4.41  through  4.44  is  worse  still  in  that  the  MFC 
controller,  and  most  likely  any  other  controller,  cannot  even  maintain  stability.  The 
stabilator  is  considered  failed  at  its  most  downward  saturation  point,  causing  the  aircraft 
to  pitch  over  and  quickly  depart  the  linear  region  represented  by  the  state  space  model. 

4.2  Pitch  Pointing 

Pitch  pointing  is  a  maneuver  in  which  the  flight  path  angle  and  the  forward 
velocity  are  held  constant  while  the  pitch  angle  tracks  a  setpoint.  For  a  fighter  aircraft, 
pitch  pointing  is  important  in  that  it  provides  “point-and-shoof  ’  capability  or,  in  other 
words,  the  ability  to  track  a  target  without  altering  the  flight  path  of  the  aircraft. 
Furthermore,  as  demonstrated  in  [4,7],  it  is  a  relatively  effortless  maneuver  that  does  not 
require  a  significant  amount  of  control  power  to  accomplish.  Similar  to  the  simulations 
of  the  vertical  translation  maneuver,  the  aircraft  is  commanded  to  begin  tracking  its 
setpoint  at  time  t  =  0  and  then  commanded  to  return  to  its  equilibrium  condition  at  t  =  5.0. 
Unless  otherwise  specified,  the  tracking  weight  is  Ry  -  diag(10‘',10'',10),  and  the  control 
weight  is  =  (lO'^)S^S.  The  sample  time  is  0.5  seconds. 
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Time  (sec) 

Figure  4.43  VT  Control  Deflections  (Failure:  SS  at  t=0) 


Time  (sec) 

Figure  4.44  VT  Control  Deflections  (Failure:  SS  at  t=0) 


The  unconstrained  pitch  pointing  maneuver  shown  in  Figures  4.45  through  4.48 
demonstrates  very  clearly  that  only  minimal  control  usage  is  necessary  to  successfully 
complete  this  maneuver;  the  only  constraint  violation  is  the  rate  constraint  on  the  trailing 
edge  flaps.  Implementing  the  constraints,  as  shown  in  Figures  4.49  through  4.52,  does 
not  significantly  affect  performance,  but  does  increase  overall  control  usage  to  achieve 
the  same  level  of  tracking.  During  the  maneuver,  the  thrust  vectoring  vanes  and  the 
stabilator  are  primarily  responsible  for  establishing  the  pitch  angle,  while  the  ailerons  and 
trailing  edge  flaps  deflect  upward  to  dissipate  the  excess  lift  caused  by  the  increased  angle 
of  attack.  Similar  to  the  vertical  translation  maneuver,  the  thrust  vectoring  also  appears  to 
play  a  role  in  minimizing  deviations  of  the  forward  velocity  from  its  equilibrium  value. 

Because  control  deflections  are  relatively  small  in  this  maneuver,  additional 
failure  scenarios  similar  to  those  shown  in  the  section  dealing  with  vertical  translation  do 
not  provide  additional  insight  into  the  problem  of  dealing  with  actuator  failures. 

However,  the  tracking  weights  chosen  for  the  pitch  pointing  maneuver  do  serve  to 
provide  a  useful  example  of  how  weightings  chosen  for  optimal  tracking  may  not  be 
ideally  suited  to  the  aircraft  following  the  failure  of  a  control  element.  For  instance,  in 
Figures  4.53  through  4.56  the  aircraft  suffers  a  simulated  multiple  failure  of  the  stabilator 
at  t=0.5  seconds  and  the  ailerons  at  t=l  .0  second.  Given  the  small  deflections  at  which 
these  failures  occur,  the  thrust  vectoring  vanes  and  the  trailing  edge  flaps  easily 
compensate  for  the  malfunctions,  and  the  initial  setpoint  tracking  performance  is  only 
slightly  worse  in  terms  of  settling  time.  Where  the  difficulty  arises  is  when  the  aircraft  is 
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1.2 


Time  (sec) 

Figure  4.49  PP  Output  Response  (Nominal  Constraints) 


Figure  4.50  PP  Output  Response  (Nominal  Constraints) 


Scaled  Inputs  Scaled  Inputs 


Figure  4.51  PP  Control  Deflections  (Nominal  Constraints) 


Figure  4.52  PP  Control  Deflections  (Nominal  Constraints) 
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1.2 


Time  (sec) 


Figure  4.53  PP  Output  Resp.  (Failures:  SS  at  t=0.5;  AS  at  t=1.0) 


Time  (sec) 


Figure  4.54  PP  Control  Defl.  (Failures:  SS  at  t=0.5;  AS  at  t=1.0) 
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Time  (sec) 

Figure  4.55  PP  Control  Defl.  (Failures:  SS  at  t=0.5;  AS  at  t=1.0) 


Time  (sec) 

Figure  4.56  PP  Control  Defl.  (Failures:  SS  at  t=0.5;  AS  at  t=1.0) 


commanded  to  return  to  its  original  setpoint.  The  response  to  this  command  is  oscillatory 
and  undesirable  from  a  handling  qualities  perspective.  The  cause  of  this  oscillatory 
behavior  can  be  linked  to  the  fact  that  the  tracking  weights  for  the  forward  velocity  and 
flight  path  angle  are  three  orders  of  magnitude  greater  than  the  tracking  weight  for  the 
pitch  angle.  Thus,  the  pitch  angle  is  given  the  least  priority  when  it  comes  to  minimizing 
its  perturbations  from  its  designated  setpoint.  To  remedy  this  situation,  we  relax  the 
weights  on  the  forward  velocity  and  flight  path  angle  and  slightly  increase  the  weight  on 
the  pitch  angle  at  t=5.0  seconds:  =  diag(20,20,4).  The  effect  on  the  output  response  is 
a  significant  decrease  in  the  magnitude  of  the  oscillations  and  only  a  slight  increase  in  the 
perturbations  of  the  flight  path  angle  and  the  forward  velocity  from  their  trim  values. 
Figures  4.57  through  4.60  contain  the  control  deflections  and  the  output  response  for  this 
case. 

4.3  Direct  Lift 

The  direct  lift  maneuver  is  one  in  which  the  angle  of  attack  is  held  constant  while 
the  flight  path  angle  and  the  pitch  angle  vary.  The  purpose  of  the  direct  lift  simulations 
presented  in  this  thesis  is  to  explore  how  varying  the  prediction,  control,  and  optimization 
horizons  affects  the  performance  of  the  controller.  In  all  cases,  the  basic  stability 
criterion  presented  in  Section  2.3.1  is  used  to  ensure  stability  of  the  resulting  closed  loop 
system.  The  control  and  tracking  weights  are  selected  to  be  =  (25,20,1)  and 
=  (1.2e-4)  for  every  simulation.  The  sample  time  is  t  =  0.5  seconds. 
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Figure  4.59  PP  Control  Deflections  (Relaxed  Weights) 
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Figure  4.60  PP  Control  Deflections  (Relaxed  Weights) 
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Displayed  in  Figures  4.61  through  4.64,  the  first  test  case  employs  the  standard 
horizon  lengths  used  throughout  this  thesis:  r  =  5,  p  =  q  -  20.  Like  the  vertical  translation 
maneuvers  shown  in  Section  4.1,  direct  lift  requires  considerable  control  power  to 
accomplish.  During  the  initial  setpoint  tracking,  the  thrust  vectoring  vanes  reach  their 
upward  saturation  point,  and  the  ailerons  reach  their  downward  saturation  point. 
Additionally,  both  the  trailing  edge  flaps  and  the  throttle  are  constrained  by  their 
respective  rate  limits. 

Increasing  the  optimization  horizon  to  r  =  10  has  no  discernible  effect  on  system 
performance.  Figures  4.65  through  4.68  indicate  that  the  control  deflections  and  the 
output  response  of  this  case  are  exactly  identical  to  the  nominal  case  of  r  =  5,  p  =  q  =  20. 
Note  that  ensuring  system  stability  with  an  optimization  horizon  increase  requires  that  the 
prediction  and  control  horizons  also  be  lengthened:  r  =  10,  p  =  q  =  25.  Similarly,  Figures 
4.69  through  4.72  show  that  an  additional  augmentation  of  the  optimization  horizon  to  r  = 
1 5  is  equally  ineffective  in  altering  system  tracking  performance. 

The  next  two  scenarios  involve  fixing  the  optimization  horizon  at  r  =  5  and  then 
varying  the  prediction  horizon  independently.  In  Figures  4.73  through  4.76  we  see  that 
increasing  the  prediction  horizon  to  p  =  25  has  no  effect  on  system  performance.  All 
control  deflections  are  identical  to  those  of  the  nominal  case,  and  the  output  response  is 
equally  unaffected.  Figures  4.77  through  4.80  show  that  a  further  increase  of  the 
prediction  horizon  to  p  =  30  has  similarly  negligible  effects  on  the  system  tracking 
performance. 
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Figure  4.67  DL  Control  Deflections  (r  =  10;  p  =  q  =  25) 


Figure  4.68  DL  Control  Deflections  (r  =  10;  p  =  q  =  25); 
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Scaled  Inputs 


Figure  4.69  DL  Output  Response  (r  =  15;  p  —  q  -  30) 
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Figure  4.70  DL  Control  Deflections  (r  =  15;  p  =  q  =  30) 
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Figure  4.78  DL  Control  Deflections  (r  =5;  p  =  30;  q  =  20) 
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Time  (sec) 

Figure  4.79  DL  Control  Deflections  (r  =  5;  p  =  30;  q  =  20) 


Time  (sec) 

Figure  4.80  DL  Control  Deflections  (r  =  5;  p  =  30;  q  =  20); 


s  o  Conclusions  and  Directions  for  Future  Research 


5.1  Conclusions 

As  shown  in  the  vertical  translation  simulations,  MPC  has  the  capability  to 
redistribute  control  power  in  order  to  mitigate  the  effects  of  single  and  multiple  control 
surface  failures  so  long  as  the  specified  failures  are  not  too  extreme.  For  minor  failures, 
the  MPC  controller  even  maintains  near-nominal  performance.  Another  interesting  result 
demonstrated  in  this  thesis  is  how  on-line  adjustment  of  the  optimization  weightings  can 
significantly  improve  system  performance,  especially  subsequent  to  a  control  element 
failure.  Demonstrated  in  the  pitch  pointing  simulations,  the  change  of  weights  following 
a  multiple  control  element  failure  significantly  reduced  pitch  oscillations  when  the 
aircraft  was  commanded  to  move  from  its  setpoint  to  its  original  equilibrium  condition. 
Thus,  the  addition  of  a  mechanism  for  adaptively  choosing  optimization  weightings 
seems  like  a  logical  addition  to  MPC  aircraft  controllers.  Lastly,  the  direct  lift 
simulations  present  a  case  for  choosing  the  MPC  horizons  as  small  as  possible  because,  in 
this  specific  instance,  increasing  them  beyond  the  minimums  for  system  stability  only 
increases  calculation  time  and  does  not  discemibly  benefit  overall  system  performance. 

5.2  Directions  for  Future  Research 

A  next  logical  step  in  the  use  of  MPC  controllers  for  control  of  aircraft  is  to 
augment  the  straight  four  state  F-18  linearized  model  with  actuator  and  sensor  dynamics 
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in  the  attempt  to  increase  the  overall  fidelity  of  the  simulation.  Furthermore,  if  possible, 
outputs  fed  into  the  estimator  should  be  those  actually  available  from  the  F-18  HARV  or 
whatever  aircraft  is  being  studied. 

In  the  arena  of  fault-tolerant  control,  the  pitch  pointing  simulations  indicate  that 
the  addition  of  an  adaptive  optimization  weight  selection  routine  could  conceivably 
improve  system  performance.  Additionally,  only  actuator  failure  scenarios  and  not  battle 
damage  scenarios  were  studied  in  this  thesis  because  in  most  cases  battle  damage  would 
have  at  least  some  effect  on  aircraft  dynamics.  Thus,  the  addition  of  a  system 
identification  routine  linked  to  the  MFC  optimizer  would  likely  enable  the  MFC 
controller  to  update  the  prediction  equations  based  on  the  updated  plant  model  and 
maintain  system  stability. 
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Appendix  Al:  F-18  HARV  Unsealed  Model 


yo(0  =  c^x^it) 


-7.5000e-2 

-2.4050e  +  l 

0 

-3.2160e+l 

-8.8250^-4 

-1.9590e-l 

9.8960e-l 

0 

-1.5000^-4 

-1.4540e-l 

-1.8770^-1 

0 

0 

0 

l.OOOOe+0 

0 

-2.2980e-2 

0 

-7.2865e-2 

3.9300e-2 

-4.1140e-2 

1.6000e-l 

-1.8000e-4 

-1.3000e-4 

-4.3000e-4 

-4.7000e-5 

-3.0000^-4 

-3.0000^-4 

-6.7000e-3 

-7.0000e-4 

-1.2000e-2 

-5.7000e-4 

7.0000e-4 

4.7000e-4 

0 

0 

0 

0 

0 

0 

1.0000^+0  0  0  0 
0  -l.OOOOe+0  0  l.OOOOe+0 
0  0  0  l.OOOOe+0 


100 


Appendix  A2:  F-18  HARV  Scaled  Model 
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Appendix  Bl:  Example  of  Controller  Setup  Script  File  for  Matlab 


WRITTEN  BY:  Derek  W.  Ebdon 

%  Script  file  for  controller  setup,  "sltl.m" 

%  VERTICAL  TRANSLATION 
%  CONTROL  SURFACES:  No  failure 

%  This  problem  will  use  the  differencing  formulation  from  equations  (2.18)  and  (2.19) 

%  to  manipulate  the  plant  into  a  form  that  accepts  control  increment  inputs. 

%  The  flight  condition  is  as  follows: 

%  Alt.:  15,000  ft 
%  Mach:  0.24 
%  Vt  :  238.7  ft/s 
%  AOA  :  25  deg 
%  PA  :  25  deg  (pitch  angle) 

%  FPA  :  0  deg  (flight  path  angle) 

%  Define  the  sample  rate  of  the  system 

Tsample  =  0.5; 

%  Perform  the  input-output  and  state  scalings  and  then  discretize  the  plant. 

%  The  plant  states  are:  Vt,  alpha,  q,  theta 
%  The  plant  outputs  are:  Vt,  gamma,  theta 

harvopl  %  Load  the  state  space  longitudinal  model  and  the  scaling  matrices. 

ascl  =  inv(Ts)*alon*Ts;  %  Perform  the  scalings 
bscl  =  inv(Ts)*blon*Ss; 
cscl  =  inv(Rs)*clon*Ts; 
dscl  =  inv(Rs)*dlon*Ss; 

[a,b,c,d]  =  c2dm(ascl,bscl,cscl,dscl,Tsample,'zoh');  %  Discretize  using  a  zero  order  hold 

%  The  F-1 8  HARV  has  6  inputs  to  the  system  : 

%  1 )  symmetric  thrust  vectoring  vane 

%  2)  symmetric  aileron  deflection 

%  3)  symmetric  stabilator  deflection 

%  4)  symmetric  leading  edge  flap 
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%  5)  symmetric  trailing  edge  flap 

%  6)  throttle 


%  Create  the  state-space  matrices  for  a  plant  that  will  accept  input 
%  increments  instead  of  absolute  control  inputs. 

A  =  [  a  zeros(4,3) ;  c*a  eye(3)  ]; 

B  =  [  b  ;  c*b  ]; 

C  =  [  zeros(3,4)  eye(3)  ]; 

D  =  zeros(3,6); 


%  Create  a  continuous  time  model  that  accepts  control  increment  inputs  for  use  in 
Simulink. 

at  =  [  ascl  zeros(4,3) ;  cscl*ascl  eye(3)  ]; 
bt  =  [  bscl ;  cscl*bscl  ]; 
ct  =  [  zeros(3,4)  eye(3)  ]; 
dt  =  zeros(3,6); 


%  Establish  the  prediction,  optimization,  and  control  horizons. 

p  =  20;  %  prediction  horizon 

q  =  20;  %  control  horizon 

r  =  5;  %  optimization  horizon 

rho  =  20; 


%  Establish  Z  matrix  weights 

Zr  =  eye(6); 

Z1  =  eye(6); 


%  Establish  R  optimization  weights 

Ry=  [5000;  0200  ;00  le4]; 

Ru  =  (le-6)*Ss'*Ss; 
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%  Determine  the  number  of  control  inputs,  outputs,  and  states. 


kappa  =  size(A,l);  %estimator  states 
xi  =  size(B,2);  %control  inputs 
eta  =  size(C,l);  %system  outputs 


%  Establish  the  setpoint  if  not  performed  in  Simulink. 
%  sinf  =  [0  1  0]'; 


%  Establish  the  constraint  matrices:  LI,  Mm,  Nn.  These  matrices  implement  the 
%  constraints  at  time  k  only. 

LI  =  [  -eye(xi)  zeros(xi,xi*(q-l)) ;  eye(xi)  zeros(xi,xi*(q-l))  ]; 

Mm  =  [  -eye(xi)  zeros(xi,xi*(q-l))  ;  eye(xi)  zeros(xi,xi*(q-l))  ]; 

Nn=[]; 


%  Establish  physical  constraint  limits:  1,  m,  n 
%  The  rate  constraints  assume  a  0.5  second  sample  time. 


1  = 

[4.0 

2.0 

1.176 

0.500 

0.300 

0.556 

4.0  2.0  1.176 

0.500 

0.300 

0.556]'; 

m  = 

n  = 

[1.0 

[]; 

1.0 

1.035 

0.222 

0.293 

1.722 

2.5  1.8  0.994 

2.178 

1.473 

0.981  ]'; 

12  = 

[4.0 

2.0 

1.176 

0.500 

0.300 

0.556 

4.0  2.0  1.176 

0.500 

0.300 

0.556]'; 

m2  = 

[1.0 

1.0 

1.035 

0.222 

0.293 

1.722 

2.5  1.8  0.994 

2.178 

1.473 

0.981  ]'; 

n2  = 

[]; 

13  = 

[4.0 

2.0 

1.176 

0.500 

0.300 

0.556 

4.0  2.0  1.176 

0.500 

0.300 

0.556]'; 

m3  = 

[1.0 

1.0 

1.035 

0.222 

0.293 

1.722 

2.5  1.8  0.994 

2.178 

1.473 

0.981  ]'; 

n3  = 

ri; 

%  Calculate  the  gains,  Fr  and  L,  that  place  the  poles  of  the  plant  and  the  estimator  at  the 
%  origin. 

[Fr,L]=gains(A,B,C,D); 
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%  Calculate  the  state  space  matrices  for  the  stabilizing  inner  loop  controller. 
[  Ax,Bx,Cx,Dx,  Ay,By,Cy,Dy,FO,F  1  ,G,H]=controller(A,B»C,Fr,L,Zl,Zr) ; 


%  Formulate  the  prediction  matrices. 

[cF,cG,cGinf,cH,cFdel,cGdel,cGdelinf,cHdel,cFu,cGu,cGuinf,cHu]= 
pmats(kappa,xi,eta,p,q,r,rho,FO,F  1  ,G,H,Fr,Zr); 


%  Calculate  the  matrices  necessary  to  perform  the  quadratic  programming  optimization. 
[S,cB,cD,cE]= 

qpconstants(C,Ry,Ru,cF,cG,cGinf,cH,cFdel,cGdel,cHdel,cGdelinf,cFu,cGu,cHu, 

cGuinf,Ll,Mm,Nn,p,q,r,rho); 


%  Calculate  the  vector  of  constraints  across  the  MFC  horizons. 
[Il,mm,nn,112,nim2,nn2,113  ,mm3  ,nn3]=constraint(l,m,n,12,m2,n2,13  ,m3  ,n3  ,p,q) ; 


%  Calculate  setpoint  trajectory  and  final  value  of  “quasi-reference”  signal,  v. 

%[vinf]  =  vinfy(A,B,C,Fr,L,Zr,sinf,r,rho); 

%[s]  =  setpoint(sinf,p); 


save  A  A 
save  B  B 
save  C  C 
save  Fr  Fr 
save  L  L 
save  Zr  Zr 
save  rho  rho 
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Appendix  B2:  Kucera  Algorithm  for  Pole  Placement 


WRITTEN  BY:  Sharon  A.  Heise 

function  [Fr,G]=gains(A,BjC,D); 

%  to  make  A+BFr  and  A+LC  stable 
%  Q=C'*Ry*C; 

%  [Fr,ricsol,evalsl]=dlqr(A,B)Q»Ru); 

%  Fr=-Fr; 

%  [L]=dlqe(A,eye(size(A)),C,eye(size(A)),eye(size(C,l))); 
%  L=-L; 

%  to  place  poles  of  A+BFr  and  A+LC  at  zero 


% 

%  just  to  remember  the  original  a,b,c,d 
% 

ra=A; 

rb=B; 

rc=C; 

rd=D; 

kappa=size(A,l); 

xi=size(B,2); 

emp=eye(xi); 

debug=0; 


X=[]; 

Y=[]; 

Yn=[]; 

XY=[]; 


if  (rank(ctrb(A,B))'^kappa) 
error('[A,B]  is  not  controllable'); 
end 


% 

%  since  a*a*a  is  more  accurate  than  a^3 

%  use  Anb=A'^i  *  B 

Anb=B; 
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for  i=l:kappa; 

% 

%  kucera's  C_i. 

% 

%disp(sprintf(’Iteration  start:  i=%d’,i)); 
XYA=[XYAnb]; 


% 

%  The  trouble  with  having  a  full  rank  test  like  this 
%  is  that  it  is  possible  to  add  _too_inany_  columns 
%  to  XY  using  D_l,  so  that  there  is  no  room  to  make 
%  D_2  etc  non-zero.  So  lets  just  skip  it  and  add  only 
%  one  column  maximum. 

% 


Yn=[]; 

testXY=XY; 

foundone=0; 

for  j=l  :size(Anb,2) 

newcol=Anb(;,j); 

if  ( rank([testXY  newcol])  ==  size([testXY  newcol],2) ) 
testXY=[testXY  newcol]; 
foimdone=l; 

Yn=[Yn  emp(:,j)] ; 

% 

%  Note  that,  in  order  to  not  make 
%  D_1  (say)  very  wide,  and  then 
%  have  to  make  D_2,3,4  etc  all  0, 

%  we  insert  a  break  here  so  that 
%  we  stop  searching  after  finding  one  rank 
%  addition 
% 

disp(sprintf('breaking  at  i=%d  j=%d',i,j)); 
break; 
end 
end 
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if  (foundone=0) 

XY 

Anb 

error(sprintf('Failed  to  find  a  D_%d',i)); 


XY=[XY  Anb*Yn] 

place(i)=size(Y,2)+l ; 
Y=[Y  Yn]; 

X=[X  Anb*Yn]; 
place  1  (i)=size(X,2); 


% 

%  update  Anb 
% 


Anb=A  *  Anb; 


end; 

place(i+ 1  )=size(  Y,2)+ 1 ; 
place  1  (i+ 1  )=size(X,2); 


if  (debug~=0) 


disp(sprintf('place=')); 

place 

disp(sprintf('place  1  =')); 
place! 

disp(sprintf('Rank  of  [XI Y1  X2Y2  ...]  -  %d',rank(X))); 


% 

%  D_i  is  given  by  Y(:,place(i):plaee(i+1)-1) 

% 
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%  C_i  is  given  by  A^(i-1)*B 
% 

%  [C_1D_1  C_2D_2  ...  C_iD_i]  is  given  by 
% 

%  X(:,l:placel(i)) 

% 


% 

%  some  test  code. 

% 

% 

dispC - '); 

for  i=l  :kappa 

% 

%  kucera's  Cl 
% 

Anb=A^(i-l)*B; 
disp(sprintf('Kuceras  C_%d=',i)) ; 

Anb 


% 

%  kucera's  D1 
% 

myd=Y  ( :  ,place(i):place(i+ 1 )- 1 ); 
disp(sprintf('Kuceras  D_%d=',i)); 
myd 


% 

%  [ClDl  C2D2  ...  CiDi] 

% 

disp(sprintf('[ClDl  C2D2  ...  CiDi]-)); 
myx=X(:,l  :placel(i)) 
disp(sprintf('with  rank=%d',rank(myx))) ; 

dispC - '); 

end 


end 
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L=zeros(xi, kappa); 


for  i=l  :kappa; 


% 

%  first  find  myx=[ClDl  C2D2  ...  CiDi] 
% 

templ=X(;,l:placel(i)); 


% 

%  all  we  have  to  do  to  get  temp2  is  copy  the  matrix  Y  up  to  D_i 
%  and  zero  out  its  first  few  columns. 

% 


temp2=Y(:,l  :place(i+l)-l); 

temp2(: ,  l:place(i)-l)  =  zeros(size(temp2( : ,  l:place(i)-l)  )); 
% 

%  Lbar=temp2  *pinv(temp  1 ) 

% 

%  we  are  really  just  solving  linear  equations  here  ;  far  better 
%  to  use  /  than  pseudo  inverse,  also  faster,  see  "help  slash" 
% 


Lbar  =  temp2/templ ; 


% 

%  slow  but  accurate  way  to  compute  (a-b*l)  ^  k 

pm=A-B*L; 

pnir=pm; 

ifi>l 
forj=l:i-l 
pmr=pmr  *  pm; 
end 
end 

L  =  L  +  Lbar  *  pmr; 
end; 

Fr=-L; 


no 


% 

%  test 
% 


eig(A-B  *  L) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%% 

% 

%  And  now  we  try  for  the  observer... 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%% 


A=A'; 

B=C’; 

kappa=size(A,l); 

xi=size(B,2); 

emp=eye(xi); 

debug=0; 


X=[]; 

Y=[]; 

Yn=[]; 

XY=[]; 


if  (rank(ctrb(A,B))~=kappa) 
error('[A,C]  is  not  observable'); 
end 


place=[]; 

placel=[]; 
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% 

%  since  a*a*a  is  more  accurate  than  a'^3 
%  use  Anb=A^i  *  B 


Anb=B; 
for  i=l  :kappa; 

% 

%  kucera's  C_i. 
% 


%disp(sprintf('Iteration  start:  i=%d',i)); 
XYA=[XY  Anb]; 


Yn=[]; 

testXY=XY; 

foundone=0; 

for  j=l  :size(Anb,2) 

newcol=Anb(:J); 

if  ( rank([testXY  newcol])  ==  size([testXY  newcol],2) ) 
testXY=[testXY  newcol]; 
foundone=l; 

Yn=[Ynemp(:,j)]; 

break; 

end 

end 

if  (foundone==0) 

error(sprintf('Failed  to  find  a  D_%d',i)); 
end 


XY=[XY  Anb*Yn]; 

place(i)=size(Y,2)+l ; 
Y=[Y  Yn]; 

X=[X  Anb*Yn]; 
placel  (i)=size(X,2); 
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% 

%  update  Anb 

% 

Anb=A  *  Anb; 
end; 

place(i+l)=size(Y,2)+l ; 
place  1  (i+ 1  )=size(X,2); 


if  (debug~=0) 


disp(sprintf('place=')) ; 
place 


disp(sprintf('placel=')); 
place 1 


disp(sprintfCRank  of  [XI Y1  X2Y2  ...]  =  %d',rank(X))); 


% 

%  D_i  is  given  by  Y(:,place(i):place(i+1)-1) 
% 

%  C_i  is  given  by  A'^(i-1)*B 
% 

%  [C_1D_1  C_2D_2  ...  C_iD_i]  is  given  by 
% 

%  X(:,l:placel(i)) 

% 


% 

%  some  test  code. 
% 

% 


dispC 
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for  i=l:kappa 


% 

%  kucera's  C 1 
% 

Anb=A^(i-l)*B; 
disp(sprintf('Kuceras  C_%d=',i)) ; 

Anb 

% 

%  kucera's  D1 
% 

myd=Y(:  ,place(i):place(i+ 1 )- 1 ); 
disp(sprintf('Kuceras  D_%d=',i)) ; 
myd 

% 

%  [ClDl  C2D2  ...  CiDi] 

% 

disp(sprintf(’[ClDl  C2D2  ...  CiDi]=')); 
myx=X( : ,  1  :place  1  (i)) 
disp(sprintf('with  rank=%d',rank(myx))); 

dispC - '); 

end 

end 

L=zeros(xi, kappa); 
for  i=l;kappa; 

% 

%  first  find  myx=[ClDl  C2D2  ...  CiDi] 

% 

templ=X(:,l  :placel(i)); 

% 

%  all  we  have  to  do  to  get  temp2  is  copy  the  matrix  Y  up  to  D 
%  and  zero  out  its  first  few  columns. 

% 
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temp2= Y( : ,  1  :place(i+ 1 )- 1 ); 

temp2(: ,  l;place(i)-l)  =  zeros(size(temp2( : ,  l:place(i)-l) )); 


% 

%  Lbar=temp2*pinv(temp  1 ) 

% 

%  we  are  really  just  solving  linear  equations  here  ;  far  better 
%  to  use  /  than  pseudo  inverse,  also  faster,  see  "help  slash" 
% 


Lbar  =  temp2/templ ; 


% 

%  slow  but  accurate  way  to  compute  (a-b*l)  ^  k 

pm=A-B*L; 

pmr=pm; 

ifi>l 
forj=l:i-l 
pmr=pmr  *  pm; 
end 
end 

L  =  L  +  Lbar  *  pmr; 
end; 


G=-L'; 

% 

%  test 
% 


eig(ra  +  G  *  rc) 
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Appendix  B3:  Inner  Feedback  Loop  Controller 


WRITTEN  BY:  Sharon  A.  Heise 

function  [Ax,Bx,Cx,Dx,Ay,By,Cy,Dy,FO,Fl,G,H]=controller(A,B,C,Fr,L,Zl,Zr); 

%  Controller 

Ax=A+L*C; 

Bx=L; 

Cx=inv(Zl)*Fr; 

Dx=zeros(size(Cx,l),size(Bx,2)); 

Ay=A+B*Fr+L*C; 

By=B*Zl; 

Cy=Fr; 

Dy=Zl; 

%  State  Predictor 

%  Ap=(A+B*Fr+L*C); 

%  Bp=[B*Zr;-L]; 

%  Cp=eye(size(Ap,l)); 

%  Dp=zeros(size(Cp,  1  ),size(Bp,2)); 

F0=(A+B*Fr+L*C); 

Fl=(A+B*Fr); 


%for  i=l:p-l; 

%  F=[F;A+B*Fr]; 
%end; 


G=B*Zr; 

H=-L; 
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Appendix  B4:  Function  to  Construct  Prediction  Matrices 


WRITTEN  BY:  Sharon  A.  Heise 
MODIFIED  BY:  Derek  W.  Ebdon 

function  [cF,cG,cGinf,cH,cFdel,cGdel,cGdelinf,cHdel,cFu,cGu,cGuinf,cHu]= 
pmats3(kappa,xi,eta,p,q,r,rho,F0,Fl,G,H,Fr,Zr); 

%  Note  this  is  for  constant  F 
%  State  prediction  equations 


o/^*****=i.***=i=***Stg^g  Equation  Prediction  Matrices********************** 


cF=[];  F 

for  i=l:p; 

cF=[cF;Fl^(i-l)*FO]; 

end; 

cG=zeros(p*kappa,r*xi);  G 
for  i=l:p; 
forj=l:r; 
ifi>=j; 

cG((i-l)*kappa+l:i*kappa,(j-l)*xi+l:j*xi)  =  Fl^(i-j)*G; 
end; 
end; 
end; 

cGinf=zeros(p*kappa,(rho-r)*xi);  G  “ 

for  i=l:p; 
for  j=l:rho-r; 
ifi-r  >=j; 

cGinf((i-l)*kappa+l  ;i*kappa,(i-l)*xi+l  :j*xi)=  Fl^(i-r-j)*G; 
end; 
end; 
end; 


cH=[];  H 

for  i=l:p; 

cH=[cH;Fl^(i-l)*H]; 

end; 


117 


0/^*  *********  *  *Control  Increment 


Prediction  Matrices******************** 


cFdel=[Fr]; 
for  i=l:q-l; 

cFdel=[cFdel;Fr*F  1  ^(i- 1  )*F0] ; 
end; 


cGdel=zeros(q*xi,r*xi); 
for  i=l:q; 
forj=l:r; 
ifi==j; 

cGdel((i-l)*xi+l:i*xi,(j-l)*xi+l:j*xi)=Zr; 
elseif  i  >  j; 

cGdel((i-l)*xi+l  :i*xi,(j-l)*xi+l  :j*xi)=Fr*Fl^(i-j-l)*G; 
end; 
end; 
end; 


cGdelinf=zeros(q*xi,(rho-r)*xi); 
for  i=l:q; 
forj=l:rho-r; 
ifi-r==j; 

cGdelinf((i-l)*xi+l:i*xi,(j-l)*xi+l:j*xi)=Zr; 
elseif  i-r>j; 

cGdelinf((i-l)*xi+l:i*xi,(j-l)*xi+l:j*xi)=Fr*Fl^(i-j-r-l)*G; 

end; 

end; 

end; 


cHdel=[]; 
for  i=l:q-l; 

cHdel=[cHdel  ;Fr*F  1  ^(i- 1  )*H] ; 
end; 

cHdel=[zeros(xi,eta);cHdel] ; 
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o^H«***********+*****  *  *  Jnput  Prediction  Matrices*  ********************** 


cFu-[];  F„ 

for  i=l:q; 

temp=zeros(size(Fr*FO)); 

forj=l:i; 

temp=temp+cFdel((j  - 1 )  *xi+ 1  :j  *xi, :); 
end; 

cFu=[cFu;temp]; 

end; 


cGu=zeros(q*xi,r*xi); 
for  i=l:q; 
forj=l:r; 
ifi==j; 

cGu((i-l)*xi+l  :i*xi,(j-l)*xi+l  :j*xi)=Zr; 
elseif  i>j; 

tenip=zeros(size(F  0)) ; 
for  k=j:i-l; 
tenip=temp+F  1  ^(k-j ) ; 
end; 

cGu((i-l)*xi+l  :i*xi,(j-l)*xi+l  :j*xi)=Fr*tenip*G+Zr; 
end; 
end; 
end; 


cGuinf=zeros(q*xi,(rho-r)*xi);  Gj 
for  i=l:q; 
for  j=l  :rho-r; 
ifi-r=j; 

cGuinf((i-l)*xi+l  ;i*xi,(j-l)*xi+l  :j*xi)=Zr; 
elseif  i-r>j; 

temp=zeros(size(F  1 )) ; 
for  k=j:i-r-l; 
temp=temp+F  1  ^(k-j ) ; 
end; 

cGuinf((i-l)*xi+l  :i*xi,(j-l)*xi+l  :j*xi)=Fr*temp*G+Zr; 
end; 
end; 
end; 
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cHu=[];  H„ 
for  i=l:q-l; 
temp=zeros(size(F  1 )); 
for  j=0:i-l; 
temp=temp+F  1  ^(j ) ; 
end; 

cHu=[cHu;Fr*temp*H]; 

end; 

cHu= [zero  s(xi,eta)  ;cHu] ; 
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Appendix  B5:  Function  to  Calculate  Quadratic  Programming  Matrices 


WRITTEN  BY:  Sharon  A.  Heise 
MODIFIED  BY:  Derek  W.  Ebdon 

function  [S,cB,cD,cE]=qpconstants(C,Ry,Ru,cF,cG,cGinf,cH,cFdel,cGdel, 
cHdel,cGdelinf,cFu,cGu,cHu,cGuinf,Ll,Mm,Nn,p,q,r,rho); 


%  Determine  the  number  of  inputs,  the  number  of  outputs,  the  number  of  estimator  states, 
%  and  the  number  of  constraints  acting  along  the  prediction  and  control  horizons. 

eta=size(C,l); 

kappa=size(C,2); 

xi=size(Ru,l); 

lambda=size(Ll,l); 

mu=size(Mm,l); 

nu=size(Nn,l); 

%  Formulate  the  weighting  matrices  and  the  matrix  with  p  output  matrices  along  the 
%  diagonal. 


for  i=l:p; 

cC((i-l)*eta+l  :i*eta,(i-l)*kappa+l  :i*kappa)=C;  C 
cQ((i-l)*eta+l  :i*eta,(i-l)*eta+l  ;i*eta)=Ry; 
end; 

for  i=l:q; 

cR((i-l)*xi+l  :i*xi,(i-l)*xi+l  :i*xi)=Ru; 
end; 


for  i=l;q; 

Iq=[Iq;eye(xi)];  / 
end; 


%  Construct  the  matrices  needed  for  the  quadratic  optimization 
%  Optimization  matrices: 

S=cG'*cC'*cQ*cC*cG+cGdel'*cR*cGdel;  S 
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cB  =  2*([cC*cF  cC*cH  cC*cGinf -eye(p*eta)]'*cQ*cC*cG  + 

[cFdel  cHdel  cGdelinf  zeros(q*xi,p*eta)]'*cR*cGdel);  T 

%  Constraint  matrices  for  the  quadratic  optimization.  Note  constraints  only  applied  to 
%  inputs; 

cD  =  [  Ll*cGdel ;  Mm*cGu  ];  D 

cE=[-Ll*cFdel  -Ll*cHdel  -Ll*cGdelinf  zeros(lambda,xi)  eye(lambda)  zeros(lambda,m_u) ; 
-Mm*cFu  -Mm*cHu  -Mm*cGuinf  -Mm*Iq  zeros(mu, lambda)  eye(mu)  ];  E 

xietkarp=[xi  eta  kappa  r  p]; 


%  Save  matrices  for  use  with  Simulink. 

save  S  S 
save  cB  cB 
save  cD  cD 
save  cE  cE 

save  xietkarp  xietkarp 
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Appendix  B6:  Function  to  Assign  Constraints 


WRITTEN  BY:  Derek  W.Ebdon 

%  This  function  is  for  establishing  constraints  when  three  cases  of  constraints  are 
%  specified.  Failures  may  be  simulated  by  assigning  a  different  set  of  constraints  over 
%  each  time  period. 

function  [Il,mm,nn,112,mm2,nn2,113  ,mm3  ,nn3  ]=constraint(l,m,n,12,m2,n2,13  ,m3  ,n3  ,p,q); 

%  Constraints  will  only  be  applied  at  time  k  for  this  thesis,  so  all  that  is  required  is  to 
%  assign  one  set  of  constraints  for  each  failure  scenario.  If  constraints  are  to  be  applied 
%  across  the  entire  control  and  prediction  horizons,  we  must  use  loops  to  assign  the 
%  selected  constraint  condition  to  every  time  step  along  the  respective  horizons. 

%  Remember  that  the  minumum  limits  are  all  in  the  first  half  of  the  constraint  vector 
%  (ll,mm,nn)  and  that  all  maximum  limits  are  in  the  second  half  because  of  the  method  in 
%  which  we  formulated  the  constrained  quadratic  optimization. 


11  =  1; 
mm  =  m; 
nn  =  n; 

112  =  12; 
mm2  =  m2; 
nn2  =  n2; 

113  =  13; 
mm3  =  m3; 
nn3  =  n3; 

%  Save  constraint  vectors  for  use  with  Simulink 

save  11 11 
save  mm  mm 
save  nn  nn 
save  112  112 
save  mm2  mm2 
save  nn2  rm2 
save  113  113 
save  mm3  nim3 
save  rm3  nn3 
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Appendix  B7: 

Function  to  Calculate  Far  Future  Value  of  Quasi-Reference  Signal 
Function  to  Establish  Setpoint  Trajectory 


WRITTEN  BY:  Sharon  A.  Heise 
Quasi-Reference  Signal  Calculation 

function  [vinf] =vinfy ( A3>C  ,F r ,L3r,sinf,r,rho) ; 

kappa=size(A,l); 

X  =  [A+B*Fr]; 

temp  =  zeros(kappa, kappa); 
for  i=l:2*kappa; 
temp  =  temp+X^(i-l); 
end; 

vinf  =  (C*temp*B*Zr)\sinf; 

for  i=l  :rho-r; 

vinfer  =  [vinfer;vinf]; 
end; 

vinf  =  vinfer; 
save  vinf  vinf 


Setpoint  Trajectory 

function  [s]=setpoint(sinf,p); 


s=[]; 


for  i=l:p; 
s=[s;sinf]; 
end; 

save  s  s 
save  sinf  sinf 
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Appendix  B8:  Simulink  Optimization  Function 


WRITTEN  BY:  Derek  W.  Ebdon 

%  This  m-file  performs  a  constrained  quadratic  optimization  to  find  the 
%  quasi-reference  signal  v(k)  that  is  input  to  the  plant. 

function  [v]  =  optimize(c); 

%  Load  the  optimization  matrices.  S2  and  cB2  are  used  if  the  weighting  matrices  are 
%  modified  to  improve  performance  following  an  actuator  failure. 

load  S 
load  cB 
load  cD 
load  cE 
%load  S2 
%load  cB2 

%  Load  the  constraint  vectors. 

load  xietkarp 
load  11 
load  mm 
load  nn 
load  112 
load  mm2 
load  nn2 
load  113 
load  mm3 
load  nn3 

%  Load  matrices  required  to  determine  vinf  (far  future  value  of  v(k)). 

load  A 
load  B 
load  C 
load  Fr 
load  L 
load  Zr 
load  rho 
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%  Determine  the  number  of  system  inputs,  system  outputs,  estimator  states,  and  the 
%  lengths  of  the  prediction  and  optimization  horizons. 

xi  =  xietkarp(l); 
eta  =  xietkarp(2); 
kappa  =  xietkarp(3); 
r  =  xietkarp(4); 
p  =  xietkarp(5); 

%  Accept  input  from  Simulink  in  order  to  formulate  optimization  problem. 

u  =  c(l:xi); 

sinf  =  c(xi+l:xi+eta); 

y  =  c(xi+eta+l  :xi+eta+eta); 

xhat  =  c(xi+eta+eta+l  :xi+eta+eta+kappa); 

t  =  c(xi+eta+eta+kappa+ 1 ); 

%  Calculate  far  future  value  of  v(k)  and  the  setpoint  trajectory. 

[vinf]  =  vinfy(A,B,C,Fr,L,Zr,sinf,r,rho); 

[s]  =  setpoint(sinf,p); 

%  Linear  term  in  the  quadratic  programming  problem.  Used  just  to  clean  up  presentation. 
%  T1  and  T2  represent  different  weights:  Ry  and  Ru  (Tl);  Ry2  and  Ru2(T2). 

T1  =  [  xhat'  y'  vinf  s'  ]*cB; 

T2  =  [  xhat'  y'  vinf  s'  ]*cB2; 

%  Calculate  the  right  hand  side  of  the  constraint  equations,  bl,  b2,  b3  all  represent 
%  different  constraint  configurations. 

bl  =  cE*[xhaf  y'  vinf  u'  11'  mm'  ]'; 
b2  =  cE*[xhaf  y'  vinf  u'  112'  mm2'  ]'; 
b3  =  cE*[xhaf  y'  vinf  u’  113'  mm3'  ]'; 

if  t  <  0.5  %  Implement  constraint  set  Ifor  all  times  less  than  the  time  specified. 
[x,lambda,how]  =  qp(2*S,Tl,cD,bl); 

elseif  t  <  1 .0  %  Implement  constraint  set  2  at  times  less  than  the  specified  time. 

[x,lambda,how]  =  qp(2*S,Tl,cD,b2); 
else  %  Implement  constraint  set  3  and  weight  set  2  at  all  other  times 
[x,lambda,how]  =  qp(2*S2,T2,cD,b3); 
end 
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%  Implement  the  first  input  move. 
v  =  x(l:xi); 
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Appendix  C:  Simulink  Diagram 


CO 
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